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KEY  TO  SYMROLS 
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m  r 

D  Period  of  time 
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F  External  force  vector 

F(  )  Functional 
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H  Average  runun,  meters 
avg 
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m 

i  (-1J 

I  Tsunami  intensity 

j  Integer 

k  Wave  number 

[K]  Coefficient  matrix 

L  Wavelength 

m  Integer 

n  Unit  normal  vector 

n  Manning's  coefficient 


n(  )  Frequency  of  occurrence 

N  Interpolation  function 

P  Pressure 

-> 

q  Velocity  vector 

{Q}  Total  load  vector 

r  Spherical  coordinate 

t  Time 

T  Wave  period 
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The  development  is  described  of  two  numerical  models  that  accurately 

simulate  the  propagation  of  tsunamis  to  nearshore  regions  and  the 

interaction  of  tsunamis  with  coastlines.   One  of  these  models  is  a 

finite  element  model  which  uses  a  telescoping  numerical  grid  to  cover 

a  section  of  the  Pacific  Ocean,  including  all  eight  islands  of  the 

Hawaiian  Islands.   The  second  model  is  a  finite  difference  scheme  that 

uses  four  rectilinear  grids  to  cover  most  of  the  west  coast  of  the 

continental  United  States.   The  finite  element  model  solves  linear  and 

dissipationless  long-wave  equations.   Such  equations  govern  nearshore 

propagation  in  the  Hawaiian  Islands  since  the  short  continental  shelf 

of  the  islands  limits  the  time  available  for  nonlinearities  and 

dissipation  to  cause  significant  effects.   The  finite  difference 

model  solves  nonlinear  long-wave  equations  that  include  bottom  stress 

terms  which  are  important  for  the  long  continental  shelf  of  the  west 


coast  of  the  United  States.   Both  models  are  verified  by  hindcasting 
actual  historical  tsunamis  and  comparing  the  numerical  model  calcu- 
lations with  tide-gage  recordings.   A  frequency  of  occurrence  analysis 
of  tsunami  elevations  at  the  shoreline  in  the  Hawaiian  Islands  is 
described.   This  analysis  is  based  upon  local  historical  data  with  the 
finite  element  model  used  to  interpolate  between  historical  data  re- 
corded during  the  period  of  accurate  survey  measurements  since  1946. 
Historical  data  recorded  at  Hilo,  Hawaii,  and  dating  to  1337  is  used 
in  conduction  with  data  recorded  since  1946  in  Hilo  and  throughout  the 
Hawaiian  Islands  to  reconstruct  elevations  at  locations  in  the  islands 
lacking  data  prior  to  1946.   Frequency  of  occurrence  curves  are  deter- 
mined for  the  entire  coastline  of  the  Hawaiian  Islands  using  these 
reconstructed  elevations.   Since  most  of  the  west  coast  of  the  United 
States  lacks  local  data  of  tsunami  activity,  a  frequency  of  occurrence 
analysis  of  tsunami  elevations  is  based  upon  historical  data  of  tsunami 
occurrence  in  tsunamigenic  regions  in  addition  to  numerical  model  cal- 
culations.  A  generation  and  deep-ocean  propagation  numerical  model  is 
used  to  propagate  tsunamis  with  varying  intensities  from  locations 
throughout  the  Aleutian-Alaskan  and  Peru-Chile  regions  to  a  water  depth 
of  500  meters  off  the  west  coast.   The  nearshore  finite  difference  model 
propagates  these  tsunamis  from  the  500  meter  depth  to  shore.   The 
frequency  of  occurrence  of  combined  tsunami  and  astronomical  tide 
elevations  is  determined  by  an  analysis  involving  the  numerical  super- 
position of  tsunamis  and  local  tides. 


CHAPTER  I 

INTRODUCTION 

Tsunamis 

Of  all  water  waves  that  occur  in  nature,  ore  of  the  most  destruc- 
tive is  the  tsunami.   The  term  "tsunami,"  originating  from  the  Japanese 
words  "tsu"  (harbor)  and  "nami"  (wave) ,  is  used  to  describe  sea  waves 
of  seismic  origin.   Tectonic  earthquakes,  i.e.,  earthquakes  that  cause 
a  deformation  of  the  sea  bed,  appear  to  be  the  principal  seismic  mech- 
anism responsible  for  the  generation  of  tsunamis.   Coastal  and  submarine 
landslides  and  volcanic  eruptions  also  have  triggered  tsunamis. 

Tsunamis  are  principally  generated  by  undersea  earthquakes  of  the 
dip-slip  type  with  magnitudes  greater  than  6.5  en  the  Richter  scale  and 
focal  depths  less  than  59  kilometers.   They  are  very  long  period  waves 
(5  minutes  to  several  hoars)  of  low  height  (a  few  feet  or  less)  when 
transversing  water  cf  oceanic  depth,   Consequently,  they  are  not  dis- 
cernible in  the  deep  ocean  and  go  unnoticed  by  ships.   Tsunamis  travel 
at  the  shallow-water  wave  celerity  equal  to  the  square  root  of  the 
acceleration  due  to  gravity  times  water  depth  even  in  the  deepest  ocean; 
because  of  their  very  long  wavelengths.   This  speed  of  propagation  can 
be  in  excess  of  500  miles  per  hour  in  the  deep  ocean. 

When  tsunamis  approach  a  coastal  region  where  the  water  depth 
decreases  rapidly,  wave  refraction,  shoaling,  and  bay  or  harbor  reso- 
nance may  result  in  significantly  increased  wave  heights.   The  great 


period  and  wavelength  of  tsunamis  preclude  their  dissipating  energy 
as  a  breaking  surf;  instead,  they  are  apt  to  appear  as  rapidly  rising 
water  levels  and  only  occasionally  as  bores. 

Over  500  tsunamis  have  been  reported  within  recorded  history. 
Virtually  all  of  these  tsunamis  have  occurred  in  the  Pacific  Basin. 
This  is  because  most  tsunamis  are  associated  with  earthquakes,  and  most 
seismic  activity  beneath  the  oceans  is  concentrated  in  the  narrow 
fault  zones  adjacent  to  the  great  oceanic  trench  systems  that  are 
predominantly  confined  to  the  Pacific  Ocean. 

The  loss  of  life  and  destruction  of  property  due  to  tsunamis  have 
been  immense.   The  Great  Iloei  Tokaido-Nanhaido  tsunami  of  Japan  killed 
30,000  people  in  1707.   In  186S,  the  Great  Peru  tsunami  caused  25,000 
deaths  and  carried  the  frigate  U.S.S.  Waterlee  1,500  feet  inland. 
The  Great  Meiji  Sanriku  tsunami  of  1896  killed  27,122  persons  in  Japan 
and  washed  away  over  10,000  houses. 

In  recent  times,  three  tsunamis  have  caused  major  destruction  in 
areas  of  the  United  States.   The  Great  Aleutian  tsunami  of  1946  killed 
173  persons  in  Hawaii,  where  runup  heights  as  great  as  55  feet  were 
recorded.  The  1960  Chilean  tsunami  killed  330  people  in  Chile,  61  in 
Hawaii,  and  199  in  distant  Japan.   The  most  recent  major  tsunami  to 
affect  the  United  States,  the  1964  Alaskan  tsunami,  killed  107  people  in 
Alaska,  4  in  Oregon,  and  11  in  Crescent  City,  California,  and  caused 
over  100  million  dollars  in  damage  on  the  west  coast  of  North  America. 


Objective  and  Scope  of  Stud}' 

The  objective  of  the  present  study  is  to  determine  tsunami  eleva- 
tion frequencies  of  occurrence  at  the  shoreline  for  the  Hawaiian  Islands 
and  the  west  coast  of  the  continental  United  States,   Such  elevation  pre- 
dictions are  required  by  the  Pacific  Ocean  Division  of  the  U.  S.  Corps 
of  Engineers  for  use  in  tsunami  flood  hazard  evaluations  for  floodplain 
management  and  the  Federal  Insurance  Administration  of  the  Department  of 
Housing  and  Urban  Development  for  flood  insurance  rate  calculations. 

The  current  trend  of  general  population  migration  to  coastal 
regions  and  the  growth  of  the  tourist  industry  makes  it  essential  to 
know  possible  tsunami  inundation  levels  even  for  coastal  regions  that 
are  currently  undeveloped.   There  arc  only  a  few  isolated  locations  in 
the  United  States  that  have  sufficient  local  historical  data  of  tsunami 
activity  to  allow  accurate  predictions  solely  based  upon  historical 
data.   This  lack  of  historical  data  of  tsunami  activity  on  most  of  the 
coastline  of  the  United  States  makes  it  necessary  to  use  numerical 
models  in  addition  to  existing  historical  data  to  predict  tsunami 
elevation  frequency  of  occurrence.   However,  numerical  models  have  not 
been  successfully  applied  in  the  past  to  the  simulation  of  actual 
tsunami  nearshore  propagation  and  interaction  with  coastlines.  There 
fore,  two  numerical  models  are  developed  in  the  present  study  to 
accurately  simulate  tsunami  nearshore  activity.   Calculations  of  these 
deterministic  numerical  models  are  combined  with  probabilistic  analyses 
based  upon  historical  data  to  predict  frequencies  of  occurrence  of 
tsunami  elevations  for  the  Hawaiian  Islands  and  the  west  coast  of  the 
United  States. 


CHAPTER  II 
TSUNAMI  INTERACTION  KITH  THE  HAWAIIAN  ISLANDS 

Introduction 

Early  studies  of  tsunami  interaction  with  islands  involved  analytic 
investigations  of  the  scattering  of  monochromatic  plane  waves  by  a 
single  island  of  simple  geometric  shape.   Oner  and  Hall  (1949)  deter- 
mined the  theoretical  diffraction  pattern  for  long-wave  scattering  by 
a  circular  cylinder  in  water  of  constant  depth  by  adopting  results  from 
acoustic  theory.   They  found  qualitative  agreement  between  this  dif- 
fraction pattern  and  smoothed  runup  heights  observed  for  the  1946 
tsunami  on  the  island  of  Kauai  in  the  Hawaiian  Islands.   Hom-ma  (1950) 
determined  the  diffraction  pattern  for  long  wave  scattering  off  a  circu- 
lar cylinder  surrounded  by  a  parabolic  bathymetry  that  extended  to  a 
prescribed  distance  offshore,  beyond  which  the  depth  was  constant. 

Hydraulic  models  have  been  used  to  study  tsunami  interaction  with 
single  islands  that  are  realistically  shaped  and  surrounded  by  variable 
bathymetry.   Van  Porn  (1970)  studied  tsunami  interaction  with  Viake 
Island  using  a  1:57000  undistorted  scale  model.   Jordaan  and  Adams 
(1968)  studied  tsunami  interaction  with  the  island  of  Oahu  in  the 
Hawaiian  Islands  using  a  1:20000  undistorted  scale  model.   They 
found  poor  agreement  between  historical  measurements  of  tsunami  runup 
and  the  hydraulic  model  data.   Scale  effects  (e.g.  viscous  effects) 
and  the  effects  of  the  arbitrary  boundaries  that  confine  the  model  are 
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problems  that  are  inherent  in  hydraulic  modeling  of  tsunamis.   The  small 
heights  of  tsunamis  relative  to  deep  ocean  depths  is  also  difficult  to 
simulate.   Jordaan  and  Adams  [1968)  modeled  waves  to  a  vertical  scale 
of  1:2000;  thus,  the  waves  had  heights  ten  times  the  normal  proportion. 
Even  with  this  distortion,  waves  typically  had  heights  in  the  model  of 
only  0.3  millimeters. 

Vastano  and  Reid  (1966)  studied  the  problem  of  the  response  of 
a  single  island  to  monochromatic  plane  waves  of  tsunami  period  using 
a  finite  difference  numerical  model.   A  transformation  of  coordinates 
allowed  a  mapping  of  the  island  shore  as  a  circle  in  the  image  plane. 
The  finite  difference  solution  employed  a  grid  which  allowed  greater 
resolution  in  the  vicinity  of  the  island  than  in  the  deep  ocean. 

Vastano  and  Bernard  (1973)  extended  the  techniques  developed  by 
Vastano  and  Reid  (1966)  to  multiple- island  systems.   However,  the 
transformation  of  coordinates  technique  allows  high  resolution  only 
in  the  vicinity  of  one  island  of  a  multiple-island  system.   Thus,  for 
the  three-island  system  of  Kauai,  Oahu,  and  Niihau  in  the  Hawaiian 
Islands,  the  two  islands  of  Oahu  and  Niihau  had  to  be  represented  by 
cylinders  with  vertical  walls  whose  cross  sections  were  truncated 
wedges.   Kauai  was  represented  by  a  circular  cylinder  with  the  sur- 
rounding bathymetry  increasing  linearly  in  depth  with  distance  radially 
from  the  island  until  a  constant  depth  was  attained.   A  single  Gaussian- 
shaped  plane  wave  was  used  as  input  to  the  model.   No  comparisons  were 
made  with  historical  tsunami  data  for  the  three  islands. 

Hwang  and  Divoky  (1975)  presented  a  simulation  of  the  1964  tsunami 
at  Hilo,  Hawaii,  using  a  finite  difference  numerical  model.   The  1964 
Alaskan  tsunami  was  generated  and  propagated  across  the  deep  ocean  using 


a  finite  difference  numerical  model  employing  a  coarse  grid.   The 
Hawaiian  Islands  were  too  small  to  be  represented  by  this  numerical  grid. 
A  waveform  was  calculated  at  a  point  in  water  5000  meters  deep  that  was 
approximately  200  kilometers  from  where  Hilo,  Hawaii,  would  have  been  if 
the  Hawaiian  Islands  could  have  been  modeled.   This  waveform  was  used  as 
input  to  a  finite  difference  grid  with  fine  grid  cells  that  covered  an 
area  approximately  22  kilometers  by  b7  kilometers  in  the  immediate 
vicinity  of  Hilo.   The  water  depth  along  the  input  boundary  of  the  fine 
grid  varied  from  approximately  250  meters  to  1500  meters.   Hwang's 
approach  neglected  the  influence  of  the  Hawaiian  Islands  and  of  the 
varying  bathymetry  and  depth  changes  between  the  deep  ocean  point  and  the 
fine  grid  input  boundary  on  the  tsunami  arriving  in  Hilo.   The  numerical 
model  calculations  for  the  1964  tsunami  were  in  poor  agreement  with 
the  Hilo  tide  gage  recording  of  this  tsunami. 

A  finite  difference  model  employing  a  grid  covering  the  eight 
major  islands  of  the  Hawaiian  Island  chain  was  used  by  Bernard  and 
Vastano  (1977)  to  study 'the  interaction  of  a  plane  Gaussian  pulse  with 
the  Hawaiian  Islands.   The  square  grid  cells  were  5.5  kilometers  on  a 
side  and  close  to  the  minimum  feasible  (because  of  computer  time  and 
cost  limitations)  size  for  a  rectangular  cell  finite  difference  grid 
covering  the  major  islands  of  Hawaii.   Yet  historical  data  indicate  that 
significant  variations  of  tsunami  elevations  occur  over  distances  much 
less  than  5.5  kilometers.   The  islands  of  Hawaii  are  relatively  small 
and  very  poorly  represented  by  a  5.5  kilometer  grid.   For  example, 
Oahu  has  a  diameter  of  only  approximately  50  kilometers  and  the  land- 
water  boundary  of  the  island  has  characteristic  direction  changes 


that  occur  over  distances  much  less  than  5.5  kilometers.   The  off- 
shore bathymetery  of  the  islands  also  varies  rapidly,  with  depth  changes 
of  more  than  1500  meters  frequently  occurring  in  distances  of  5.5  kil- 
ometers.  Furthermore,  if  a  resolution  of  eight  grid  cells  per  tsunami 
wavelength  is  maintained  for  tsunami  periods  as  low  as  15  minutes,  a 
5.5  kilometer  grid  cannot  be  used  to  propagate  waves  into  depths  less 
than  approximately  150  meters.   The  processes  that  cause  significant 
wave  modifications  and  subsequent  rapid  variations  of  elevations  along 
the  coastline  (that  are  known  from  historical  observations  to  occur 
during  tsunami  activity  in  the  Hawaiian  Islands)  probably  occur  within 
this  region  extending  from  water  at  a  depth  of  150  meters  to  the 
shoreline. 

In  the  present  study,  a  finite  element  numerical  model  is  used  to 
propagate  tsunamis  from  the  deep  ocean  to  the  shoreline  of  the  Hawaiian 
Islands.   Since  finite  element  techniques  allow  dramatic  changes  in 
element  sizes  and  shapes,  this  model  has  great  advantages  relative  to 
finite  difference  numerical  models  in  accurately  representing  land 
shapes,  ocean  bathymetry,  and  tsunami  waveforms.   Element  size  can  be 
large  in  deep  water  where  bathymetric  variations  are  gradual  and  wave 
lengths  are  long.   As  a  wave  enters  shallow  water  its  length  decreases 
and  the  elements  of  the  grid  can  be  telescoped  to  smaller  sizes  with  no 
accompanying  less  of  resolution. 

Governing  Equations 

Typical  tsunamis  have  lengths  much  greater  than  the  water  depths 
over  which  they  propagate;  for  such  waves,  fluid  motions  are 


approximately  two-dimensional,  i.e.,  vertical  fluid  velocities  and 
accelerations  are  small  in  comparison  with  horizontal  velocities  and 
accelerations,  respectively. 

Consider  a  small  element  of  fluid  as  shown  in  Figure  1.   The  rate 
of  increase  of  volume  of  the  element  is  given  by 


—  {(n+h)  dxdy}  =  j-  dxdy 


where, 

n  =  a  free-surface  elevation, 

h  =  the  water  depth, 

t  =  time, 
and  (x,y)  are  Cartesian  coordinates. 


(2.1) 


r,4h 


Figure  1.   Element  of  Fluid 

In  Equation  (2.1),  it  is  assumed  that  h=h(x,y),  i.e.   h  has  no  time 
dependence.   The  net  rate  of  volume  flux  is  given  in  difference  form  by 


{u(n+h)>   A  ,  dy  -  {u(n+h)}  dy 

+  {v(n+h)>   ,   dx  -  {v(n+h)>  dx 
y+dy  )' 


=  dxdy  {-7—  u(n+h)  +  3y 

oX 

=  dxdy  V  •  q(n+h) 


v(n+h)l  +  0(dx2dy  ) 


(2.2) 


where  u,v  are  horizontal  velocity  components,  w  is  the  vertical 
velocity  component  (w=0  here),  and  q  =  (u,v,w)  is  a  velocity  vector 

(a  Taylor  expansion  is  used  in  liquation  (2.2)  and  terms  of  order 

2  2 
dx  dy   have  been  dropped).   Equating  the  rate  of  increase  of  volume  to 

the  rate  of  volume  flux  yields  the  following  continuity  equation: 


f£  +  V  •  q(mh)  =0  (2.3) 

at 

The  governing  momentum  equation  for  incompressible  and  constant 
viscosity  fluid  flow  is  the  Navier-Stokes  equation  in  the  form 

p  ||-  +  q  •  Vq  =  pF  -  Vp  +  yV2q  (2.4) 

where 

p  =  the  fluid  pressure, 

p   =  the  mass  density, 

u  -   the  dynamic  viscosity  constant, 

F  =  an  external  force  vector. 

Since  vertical  fluid  velocities  and  accelerations  are  small  and 
assuming  that  gravity  is  the  only  external  force  and  the  fluid  is  in- 
viscid  so  that  the  last  term  of  Equation  (2.4)  can  be  dropped,  Equation 
(2.4)  becomes  in  component  form 


3u      3u      3u     1  3p  ,.  .. 

—  +  u  —  +  v  —  = -^-  (2.5) 

3t      3x      3y     p  9x 


3v      3v      3v     1  3p  ,-,,-, 

-—     +   u  t-  +  v  —  = -^  (2.6) 

3t      3x      3y     P  3y 
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When  integrated,  Equation  (2.7)  yields  (with  atmospheric  pressure 
equal  to  zero) 

p  =  -pg(S-Ti)  (2.8) 

i.e.,  the  pressure  is  hydrostatic,  and 

Vp  =  pgVn  (2-9) 

Equations  (2.5)  and  (2.6)  can  be  recombined  to  yield 

||+  (q   '  V)q  =    -gVn  "(2.10) 

Since  along  most  of  the  path  of  propagation  of  a  tsunami  the 
wave  amplitude  is  much  less  than  the.  water  depth,  Equation  (2.3)  can  be 
approximated  as 

|£+  V  •  (hq)  =  0  (2.11) 

It  can  also  be  shown  that  the  second  term  in  Equation  (2.10)  is 
generally  less  than  the  other  terms  by  a  factor  of  the  order  of  n/h  . 
Therefore,  for  n/h  snail,  Equation  (2.10)  is  approximated  by 


8q 

9t 

Taking  a  time  derivative  of  Equation  (2.12)  and  substituting  Equation 


g7n  (2.12) 


(2.11)  into  the  right-hand  side  of  the  subsequent  equation  yields 

i-|  =  gV(V  •  (hq))  (2.13) 

dt 

Equation  (2.13)  is  known  as  the  linear  long  wave  equation.   It  is  a 

governing  equation  for  small  amplitude  long  waves  in  an  inviscid  and 

incompressible  variable  depth  fluid  and  is  valid  when  the  assumptions 

made  in  its  derivation  remain  valid. 

Assuming  that  the  fluid  is  irrotational   q  can  be  expressed 

as  a  gradient  of  a  scalar  function  <I>(x,y,t)  ,  known  as  the 
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velocity  potential,  i.e, 

q  =  VC  (2.14) 

Assuming  all  motions  to  be  simple  harmonic  in  time  yields 

Hx,y,t)  =  Re  [OC-x.yJe-1^]  (2.15) 

where  Re  [  ]   indicates  the  real  part  of  the  braced  quantity,  and  w 

is  the  frequency  of  the  motion.   Substituting  Equations  (2.1.4)  and 

(2.15)  into  (2.13)  yields 

2 
7  •  (hVO)  +  ~   i  -  0  (2.16) 

Equation  (2.16)  is  the  well  known  generalized  Helmholtz  equation. 

The  wave  amplitude   r,  can  be  related  to  the  velocity  potential   $ 

by  evaluating  the  Bernoulli  equation  at  the  free  surface.   The  Bernoulli 

equation  is  given  by 

3 ;   P       1 1   1 2 

_  =  _  -  g7  +    vol   =  0  ,  (2.17) 

at    p        2 '   '  v    ' 

however;  the  last  term  in  Equation  (2.17)  is  generally  smaller  than  the 
others  by  a  factor  of  the  order  of  n/h  and  can  be  neglected.   Also, 
taking  p=0  at  the  free  surface  where  I-n  ,  Bernoulli's  equation 
becomes  to  the  same  order  of  approximation  as  Equation  (2.16) 

Several  assumptions  are  made  in  the  derivation  of  Equations  (2.16) 
and  (2.18)  and  the  adequacy  of  these  assumptions  for  the  particular 
problem  of  tsunami  interaction  with  the  Hawaiian  Islands  must  be  con- 
sidered.  An  important  factor  in  this  consideration  is  the  fact  that  the 
continental  shelf  of  the  Hawaiian  Islands  is  exceptionally  short. 
Thus,  during  most  of  the  time  of  tsunami  propagation,  fluid  velocities 
are  much  less  than  wave  propagation  velocities  since  wave  amplitudes 
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arc  much  loss  than  water  depths  and  bottom  stress  forces  as  well  as 
nonlinearities  have  little  time  to  cause  significant  effects. 

Ursell  (1953)  and  Hammack  (1972)  have  shown  that  nonlinear  effects 
grow  with  time.   However,  for  a  typical  large  tsunami,  such  as  the  1964 
Alaskan  tsunami}  Hammack  and  Segur  (1977)  have  shown  that  neither  non- 
linearity  nor  frequency  dispersion  have  any  significant  effect  on  the 
lead  wave  as  it  propagates  across  the  ocean.   In  addition,  when  the  lead 
wave  reaches  a  continental  shelf  a  certain  time  must  pass  before  non- 
linearity  or  frequency  dispersion  become  significant.   For  a  typical 
large  tsunami,  Hammack  and  Segur  (1977)  show  that  linear  and  nondis- 

persive  theory  holds  for  the  lead  wave  over  a  distance  of  approximately 
2 

±L^—  (a  is  a  wave  amplitude).  Taking  the  average  depth  of  the  region 
from  the  edge  of  the  continental  shelf  (600  foot  depth)  to  shore  as  300 

feet  and  the  average  amplitude  as  3  feet  (the  amplitude  becomes  large 

?h2 
only  near  shorej  yields  =  12  miles.   This  distance  is  much  greater 

than  the  length  of  the  continental  shelf  of  the  Hawaiian  Islands 

(usually  less  than  1  mile).   Since  the  lead  wave  is  followed  by  waves 

which  have  been  reflected  by  land  areas  in  the  generation  region,  linear 

and  nondispersive  theory  should  be  adequate,  at  least  for  the  major 

waves  of  a  tsunami. 

Finite  Element  Model 


Figure  2  shows  the  problem  of  tsunami  interaction  with  islands 
with  the  process  divided  into  two  domains.   One  domain  (A)  includes  the 
islands  and  the  surrounding  variable  depth  region.   The  second  domain 
(B)  is  a  constant  depth  region  extending  to  infinity.   This  constant 
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Figure  2.   Regions  of  Computation 


14 


depth  region  is  realistic  since  the  Pacific  Ocean  has  a  remarkably 
constant  depth  beyond  the  immediate  region  of  the  Hawaiian  Islands. 
The  boundary  value  problem  expressed  in  terms  of  these  domains  becomes 

the  following: 

2 

V  .,  (hr:)  +  —  5  =  0       in  Region  A  (2.19) 

2 

V  •   (Vo)  +  !L-  <?  =  0       in  Region  B  (2.20) 

gh 

$A  =  0B  on  3A  (2.21) 

A         A 

~  =   0  on  9C  =  3C  .+3CL 3C .    (2.23) 

3n  12      J 

lim  (r)1/2  C|pik)*s  =  0  (2.24) 

r  ->   co 

where 

j  =  number  of  land  masses, 

r  =  a  spherical  coordinate, 

k  =  the  wave  number, 
§  =  the  scattered  wave  velocity  potential, 

i  =  (-D    • 

Equations  (2.21)  and  (2.22)  express  the  continuity  of  the  velocity 
potential  and  its  derivative  along  the  boundary  separating  the  two 
domains  and  Equation  (2.23)  expresses  the  impermeability  of  solid  bound- 
aries.  Equation  (2.24)  is  the  Sommerfeld  radiation  condition  that 
requires  the  scattered  wave  to  be  an  outgoing  wave  at  infinity. 

A  calculus  of  variations  approach  can  be  used  to  obtain  an  Euler- 
Lagrange  formulation  of  boundary  value  problems.   This  variational 
formulation  is  based  upon  the  principle  that  of  all  possible  displace- 
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merit  configurations  a  body  can  assume  that  satisfy  compatibility  and 
the  constraints  of  kinematic  boundary  conditions,  the  configuration 
satisfying  equilibrium  makes  the  potential  energy  a  minimum.   The  Euler- 
Lagrange  formulation  of  boundary  value  problems  has  often  been  used  in 
classical  dynamics  (Hamilton's  principle  of  least  action)  and  solid 
mechanics  (structural,  solid,  and  rock  mechanics).   Of  course,  it  is  not 
possible  to  obtain  analytical  solutions  for  many  engineering  problems. 
The  finite  element  method  is  a  numerical  method  first  introduced  by 
Turner  et  al.  (1956)  (although  some  of  the  underlying  ideas  were  dis- 
cussed by  Courant  (1943))  which  has  found  wide  application  in  solid 
mechanics  for  problems  based  upon  an  Euler- Lagrange  formulation  involving 
complex  material  properties  or  boundary  conditions  (see  Zienkiewicz, 
1971  and  Desai  and  Abel,  1972).   The  finite  element  method  is  a  discrete 
approximation  procedure  applicable  whenever  a  variational  principle  can 
be  formed.   This  method  has  a  very  short  history  in  fluid  mechanic 
applications  despite  advantages  over  conventional  numerical  schemes  such 
as  finite  difference  methods.   The  first  application  of  the  finite 
element  method  to  fluid  mechanic  problems  involved  a  study  by 
Zienkiewicz  and  Cheung  (1966)  of  seepage  through  porous  media. 

Chen  and  Mei  (1974)  studied  the  problem  of  the  forced  oscillation 
of  a  small  basin  protected  by  a  breakwater  and  containing  floating 
nuclear  power  plants  by  using  the  finite  element  method.   They  dis- 
cretized  a  small  constant  depth  region  using  approximately  250  triangular 
elements  that  were  all  about  the  same  size  and  shape.   Since  small  ampli- 
tude long  waves  were  used  as  a  forcing  function,  Equation  (2.20)  was  the 
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governing  equation  (the  water  depth  being  constant  in  their  application) . 
In  the  present  study  an  approach  similar  to  that  of  Chen  and  Mei  is 
used  to  solve  the  problem  of  tsunami  interaction  with  islands  in  a 
variable  depth  region. 

For  the  finite  element  method  to  be  useful  for  tsunami  problems,  it 
is  necessary  to  take  advantage  of  two  important  properties  of  'this 
algorithm.   First,  the  finite  element  grid  which  discretizes  a  domain 
can  have  grid  cells  (elements)  of  arbitrary  size  and  shape.   Thus,  the 
grid  can  telescope  from  a  large  cell  size  in  one  section  of  a  variable 
depth  domain  to  a  very  fine  cell  size  in  another.   Secondly,  if  the 
form  of  the  spatial  variation  of  a  physical  parameter  is  known  a  priori, 
the  finite  element  method  can  employ  a  very  small  number  of  elements  in 
the  discretization  used  to  estimate  the  exact  variation.   For  example,  a 
single  finite  element  can  be  used  to  solve  the  problem  of  the  seiching 
of  a  lake  for  a  constant  depth  rectangular  lake  if  the  basic  form 
(unknown  coefficients)  of  the  solution  is  known.   The  finite  element 
method  degenerates  into  the  well-known  Rayleigh-Ritz  variational  method 
in  this  case.   This  property  of  the  finite  element  method  will  be  im- 
portant in  modeling  the  infinite  region  B  shown  in  Figure  2. 

The  variational  principle  for  the  boundary  value  problem  given 
by  Equations  (2.19)  through  (2.24)  requires  a  certain  functional 
F(-}>)  to  be  stationary  with  respect  to  arbitrary  first  variation  of  0  . 
The  first  variation  of  the  functional  F(4')  for  this  problem  is  well- 
known  and  can  be  expressed  as 
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2  3<54> 

;f  co)  =  -  //[v  •  (iwo)  ♦  £-  *]  6?  +  ;  h(iB-;AJ  H- 

A  °  "-A  A 

The  functional  can  be  derived  from  this  expression  for  its  first  vari- 
ation (Chen  and  Mei,  1974)  and  is 

?  3($  -$T) 

Ft.)  -  //  i[M«)z  - 1-  *2i  *  f  *  «VV  thH- 

A  fe        -A  A 

3(VV        3i1 
SA  A   3nA     3A  A  3nA 

3($R-$T) 
+  /  ho    gg   T    ,  (2.26) 

3A  l        SnA 
where  $   is  the  velocity  potential  of  the  incident  waves  (forcing 

function) . 

The  finite  element  method  (unlike  the.  finite  difference  method) 
discretizes  the  domain  itself  rather  than  the  governing  equations. 
The  region  under  consideration  is  partitioned  into  small  regions  known 
as  finite  elements  and  the  assemblage  of  all  such  elements  represents 
the  original  region.   Instead  of  solving  the  problem  for  the  entire  body 
in  one  operation,  the  solutions  are  formulated  for  each  constituent 
element  and  combined  to  obtain  the  solution  for  the  original  region. 
Thus,  for  example,  the  first  term  of  Equation  [2.26)  is  evaluated  on  the 
element  level  and  the  results  for  all  elements  are  summed.   This  sum  is 
equivalent  to  the  original  integral  over  the  complete  region  A. 

In  order  to  evaluate  the  first  term  in  Equation  (2.26)  on  the 
element  level,  it  is  necessary  to  approximate  the  field  variable  0 


within  each  element  by  a  polynomial  N  ,  called  a  trial  or  interpolation 
function,  in  the  coordinates   x  and  y  . 
For  each  element 

$  =  {N}  {e-  (2.27) 


where 


{£} 


I 

0 

i 

_\ 

(2.28) 


is  a  column  vector  with   0.  ,  §..  0,   the  unknov.n  values  of  the  field 

i'   3'   k 

variable  at  the  nodal  points  (vertices)  of  a  triangular  element  and 


is  a  row  vector  with 


(N)=  [N.,N  Nk] 


(a.  +  b.x  +  c.y)/2£(i=l,2,3) 

i    i     i 


a  .  -  x.  y .  -  x  .  y, 
i    k7j    y\ 


C  .  =  X  .  -  X, 

i    j    k 


(2.29) 

(2.30) 
(2.31) 
(2.32) 
(2.35) 


and  A   is  the  area  of  the  element.   The  coordinates  of  the  nodal  point 

i  are  (x.,y.)   and  the  N.   and  X,   are  similarly  expressed  by  the 

v  l  J  i  J  j        k  r        J 

cyclical  permutation  of   i,  j,  and  k  . 

A  linear  interpolation  function  can  be  used  for  tsunamis  propa- 
gating within  region  A  since  elements  within  this  region  will  always  be 
much  smaller  than  local  wavelengths  (long  waves  are  adequately  repre- 
sented as  linear  variations  over  such  short  distances).   The  convergence 
of  the  solution  is  assured  as  the  elements  are  made  infinitesimal 
(Strang  and  Fix,  1973).   Since  each  element  has  its  own  polynomial 
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expansion  independent  of  all  other  elements,  the  treatment  of  the  entire 
domain  is  systematically  handled  by  summing  the  contributions  from  each 
element. 

The  expression  $R-*T  in  Equation  (2. 26 j  is  the  scattered  wave 
velocity  potential  $„  .  Since  region  B  has  a  constant  depth,  £>  can 
be  solved  analytically.  The  governing  equation  for  this  case  is  the 
Helmholtz  equation  (Equation  (2.20)).  The  boundary  condition  at  in- 
finity requires  the  scattered  wave  obey  the  Sommerfeld  radiation 
condition  given  by  Equation  (2.24).  The  well-known  solution  to  this 
boundary  value  problem  is  as  follows: 

$  -   I  H  (kr)  (a  cos  in?  +  s  sin  m9)  (2.34) 

s    „  m       m         m 
m=0 

where 

H   are  Hankel  functions  of  the  first  kind  of  order  m, 
m 

a       and  6   are  constant  coefficients, 
m       in 

and   9   is  a  spherical  coordinate. 

In  the  finite  element  formulation,  the  infinite  region  B  can 
be  considered  to  be  a  single  element  with  the  interpolation  function 
given  by  Equation  (2.54).   Thus,  within  region  A  the  interpolation 
functions  are  linear,  but  within  region  B  the  interpolation  function 
is  an  infinite  series  involving  unknown  coefficients.   The  terms  in 
Equation  (2.54)  are  oscillatory  but  have  a  monotonically  decreasing 
modulation;  hence,  the  series  can  be  truncated  after  a  finite  number  of 
terms  depending  upon  the  desired  accuracy  (see  next  section) . 

The  required  stationarity  of  the  functional   F(e)  with  respect  to 
arbitrary  variations  of  $  at  each  nodal  point  yields  a  set  of  linear 
algebraic  equations  when  expressed  on  the  element  level  that  can  be 
solved  using  matrix  methods.   For  each  specified  incident  wave  period, 


the  velocity  potential   0  can  be  determined  at  all  nodal  points.  The 
linearized  dynamic  free  surface  boundary  condition  (Equation  (2.18))  can 
then  be  used  to  calculate  the  surface  elevation   n  . 

The  set  of  linear  algebraic  equations  evolving  from  the  element 
level  can  be  written  in  the  simple  matrix  forn 

[K]  W   =  {Q}  (2.35) 

where   [K]   is  the  large  coefficient  matrix,  {0}   is  a  column  vector 
representing  a  combination  of  nodal  unknowns  for  the  velocity  potential 
and  the  coefficient  unknowns  from  the  single  element  covering  region  B, 
and  {Q}  is  the  total  load  column  vector  related  to  the  forcing  function. 
in  the  finite  element  solution  of  structural  mechanics  problems,  one 
obtains  a  similar  equilibrium  relation  with  [K]  being  a  stiffness  matrix, 
{^}   a  nodal  displacement  vector,  and  {Q}  a  nodal  force  vector  (Desai 
and  Abel,  1972).   The  coefficient  matrix  [K]  is  symmetric,  banded,  and 
complex  (complex  notation  was  used  to  represent  waveforms,  e.g.,  the 
Hankel  functions  are  complex) .   For  the  case  of  tsunami  interaction  with 
the  Hawaiian  Islands,  [K]  will  be  exceptionally  large  and  sparse  (many 
zero  terms) .   Of  course,  only  half  of  the  banded  part  of  the  matrix 
alone  needs  to  be  manipulated  in  the  solution  process.   The  banded  part 
of  the  matrix  [K]  is  solved  using  standard  Gaussian  elimination  methods. 
The  matrix  solution  used  by  Chen  and  Mei  (1974)  was  modified  in  this 
study  to  take  advantage  of  the  sparseness  of  the  matrix  in  order  to  re- 
duce computational  time.   Also,  the  elimination  method  was  modified  so 
that  calculations  involved  only  small  blocks  of  terms  at  a  time  with  the 
remainder  of  the  matrix  kept  in  pe7%ipheral  storage.   This  matrix  parti- 
tioning was  necessary  for  the  finite  element  grid  shown  in  the  next 


21 


section  since  the  banded  part  of  the  matrix  had  over  S00,000  terms. 
Partitioning  is  possible  because  the  symmetric  coefficient  matrix  is 
positive  definite;  hence,  a  solution  is  possible  by  elimination  methods 
without  pivoting.   Without  pivoting,  elimination  performed  using  one 
row,  affects  only  the  triangle  of  elements  within  the  band  below  that 
row.   Thus,  only  two  triangular  submatrices  needed  to  be  stored  in 
computer  memory  at  a  given  time. 

Numerical  Grid 

The  finite  element  grid  covering  the  variable  depth  region  A 
is  shown  in  Figure  3.   Elements  of  the  grid  telescope  from  large  sizes 
in  the  deep  ocean  to  triangles  with  areas  as  small  as  0.5  square 
kilometers  in  shallow  coastal  waters.   The  grid  covers  an  area  of 
approximately  450,000  square  kilometers.   In  general,  finite  element 
techniques  allow  elements  to  be  any  arbitrary  shape  (e.g.  quadrilateral) 
provided  that  the  ratio  of  the  lengths  of  the  shortest  and  longest 
sides  of  an  element  is  not  extreme.   Triangular  shapes  were  used 
for  convenience  in  the  present  study. 

The  geometric  shapes  of  the  eight  islands  comprising  the  Hawaiian 
Island  chain  and  the  rapid  bathymetric  variations  surrounding  the.  islands 
are  modeled  very  precisely  by  the  finite  element  grid  of  Figure  3. 
Also,  the  number  of  node  points  along  the  shorelines  of  the  islands 
is  very  dense.   Since  wave  heights  are  calculated  at  node  points  by 
the  finite  element  model,  the  model  can  adequately  represent  the  rapid 
wave-height  variations  along  coastlines  that  are  known  to  occur 
during  tsunami  activity  in  the  Hawaiian  Islands. 


Figure   3.       Finite    I; lenient    Grid. 


The  finite  element  grid  was  drawn  by  hand  with  the  element  sizes 
dictated  by  the  local  water  depth.   The  lengths  of  element  sides  were 
allowed  to  be  no  greater  than  one-eighth  of  the  local  wavelength  for  the 
shortest  period  waves  considered.   Such  a  cell  size  results  in  a  wave 
resolution  error  of  no  more  than  one  percent  f.Mei  197S)  .   There  are 
computer  programs  (Cole  and  Riech,  1976)  that  generate  finite  element 
grids  and  number  the  nodes  to  minimize  bandwidths;  however,  such  pro- 
grams only  can  be  applied  to  fairly  homogeneous  domains.   The  rapid 
depth  variations  in  region  A  require  a  rapidly  telescoping  grid  that  is 
too  complicated  for  automatic  grid  generation  techniques. 

Each  node  of  the  finite  element  grid  must  be  numbered  and  this 
numbering  is  extremely  important  since  it  determines  the  bandwidth  of 
the  large  coefficient  matrix  discussed  in  the  last  section.   The  band- 
width of  this  matrix  is  controlled  by  the  largest  difference  (among  all 
elements)  in  the  nodal  numbers  between  two  nodes  of  the  same  element. 
The  computational  time  of  the  numerical  model  is  proportional  to  the 
cube  of  the  bandwidth  of  the  coefficient  matrix. 

For  a  simple  region  covered  by  homogeneous  elements,  the  bandwidth 
is  minimized  by  numbering  nodal  points  in  cyclic  sequence  following  a 
spiraling  pattern  '..here  one  spiral  is  shielded  by  the  following  spiral 
from  a  third  spiral.   The  grid  pictured  in  Figure  3,  however,  is  so 
large  and  complicated  that  considerable  judgement  and  effort  is  involved 
in  determining  a  numbering  pattern  that  minimizes  the  bandwidth,  thus 
resulting  in  reasonable  computational  times  for  the  numerical  model. 

The  numerical  model  also  requires  the  spatial  coordinates  of  all 
node  points  and  information  concerning  how  the  elements  are  connected 


relative  to  each  other  (connectiveness) .   The  coordinates  were  deter- 
mined using  an  electronic  digitizer  and  a  simple  computer  program 
developed  to  relate  the  internal  coordinate  system  of  the  digitizer 
to  a  prototype  coordinate  system  through  coordinate  rotation  and 
scaling.   The  elements  are  numbered  in  an  arbitrary  fashion  and  then 
connectiveness  information  assembled  for  the  numerical  model  by  listing 
the  three  nodes  of  each  element. 

The  nodal  numbering,  digitizing,  and  connectiveness  assembly 
necessarily  introduce  errors.   For  example,  the  connectiveness  assembly 
involves  determining  over  17,000  numbers  and  placing  them  on  computer 
cards.   Many  connectiveness  mistakes  do  not  result  in  error  conditions 
which  terminate  the  numerical  model's  computations.   Instead,  incorrect 
answers  are  produced  by  the  model.   To  eliminate  all  errors  in  the 
input  data,  computer  plotting  programs  were  developed  to  draw  the 
grid  and  number  nodes  and  elements  based  upon  the  connectiveness  and 
node  coordinate  information  used  as  input  to  the  numerical  model.  Such 
plots  make  errors  readily  apparent. 

The  number  of  terms  retained  in  the  truncation  of  Equation  (2.34) 
was  determined  by  trial  and  error.  The  final  number  of  terms  was  such 
that  addition  of  further  terms  produced  negligible  effects.   It  was 
found  that  this  occurs  for  terms  approximately  10  "  times  smaller  than 
the  initial  terms. 

The  computational  time  requirements  of  the  finite  element  model  for 
the  grid  shown  in  Figure  3  are  very  modest,  making  it  economically 
feasible  to  determine  the  interaction  of  arbitrary  tsunamis  with  the 
Hawaiian  Islands  (discussed  in  next  section) .   The  reason  that  the 


computationtal  time  required  by  the  finite  element  model  is  reasonably 
small  is  that  the  grid  uses  small  elements  only  in  areas  where  they 
are  necessary.   The  grid  of  Figure  3  has  approximately  2500  points, 
whereas  the  finite  difference  grid  of  Bernard  and  Vastano  (1977)  had 
26,000  grid  points.   Even  so,  some  of  the  elements  of  Figure  3  are  as 
much  as  sixty  times  smaller  than  the  finite  difference  grid  cells  used 
by  Bernard  and  Vastano. 

Model  Verification 

The  finite  element  model  is  verified  in  the  present  study  by 
comparing  numerical  simulations  of  the  1960  Chilean  and  1964  Alaskan 
tsunamis  with  tide  gage  recordings  of  these  tsunamis  in  the  Hawaiian 
Islands.   These  two  tsunamis  are  the  only  major  tsunamis  for  which 
reliable  information  exists  concerning  characteristics  of  the  ground 
uplift  which  generated  the  waves  (much  more  information  exists  for  the 
Alaskan  source) . 

A  deep  ocean  recording  of  a  tsunami  near  the  Hawaiian  Islands 
has  never  been  made.   Therefore,  prototype  wave  records  in  the  deep 
ocean  cf  the  1960  Chilean  and  1964  Alaskan  tsunamis  are  not  available 
for  use  as  input  to  verify  the  finite  element  model.   However,  a  finite 
difference  numerical  model  which  solves  the  linear  long  wave  equations 
employing  a  spherical  coordinate  grid  has  been  developed  (Hwang  et  al., 
1972),  and  used  in  several  studies  (Houston  and  Garcia,  1974;  Houston 
and  Garcia,  1975;  and  Houston  et  al . ,  1975b)  to  generate  tsunamis  and 
propagate  them  across  the  deep  ocean.   This  deep  ocean  numerical  model 
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is  used  in  this  study  to  determine  deep  ocean  waveforms  of  the  1960  and 
1964  tsunamis  that  were  needed  to  verify  the  finite  element  model. 

The  deep-ocean  numerical  model  simulates  the  uplift  deformation  of 
the  ocean  water  surface  caused  by  the  permanent  vertical  displacement  of 
the  ocean  bottom  during  an  earthquake  and  the  subsequent  propagation  of 
the  resulting  tsunami  across  the  deep  ocean.   The  permanent  deformation 
(permanent  in  the  sense  that  the  time  scale  associated  with  it  is  much 
longer  than  the  period  of  the  tsunami)  and  not  the  transient  movements 
within  the  time-history  of  the  ground  motion  is  considered  to  be  the 
important  parameter  governing  far- field  wave  characteristics.   Tran- 
sient movements  occur  for  a  period  of  the  order  of  tens  of  seconds, 
whereas,  tsunami  wave  periods  are  of  the  order  of  tens  of  minutes. 
Experimental  investigations  of  tsunami  generation  by  Ilammack  (1972)  show 
that  the  transient  ground  movements  do  not  influence  far- field  char- 
acteristics of  resulting  tsunamis  for  spatially  large  ground  displace- 
ments occurring  over  a  short  period  of  time. 

Hwang  et  al ,  (1972)  verified  the  deep  ocean  model  by  simulating 
the  1964  Alaskan  tsunami.   The  permanent  ground  deformation  in  the 
source  region  was  taken  from  Plafker  (1969)  and  a  1/4°  by  1/6°  grid  was 
used  for  the  simulation.   Since  this  grid  was  much  too  coarse  to  allow 
accurate  modeling  of  shallow-water  propagation,  a  comparison  was  made 
between  the  numerical  model  calculations  and  a  gage  recording  of  the 
1964  tsunami  at  Wake  Island  where  Van  Dorn  (1970)  had  a  pressure  trans- 
ducer in  relatively  deep  water  (800  foot  depth)  some  distance  offshore. 
The  amplitude  and  length  of  the  first  wave  of  this  tsunami  were  shown 
to  be  in  good  agreement.   (Actually,  the  grid  used  by  Hwang  et  al. 
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(1972),  did  not  extend  as  far  as  Wake  Island.   The  deep  ocean  model 
propagated  the  tsunami  to  a  point  in  the  Pacific  Ocean  and  then  simple 
refraction  techniques  were  used  to  estimate  the  change  in  the  waveform 
during  propagation  from  this  point  to  Wake  Island.) 

The  simulation  of  the  1960  and  1964  tsunamis  performed  in  this 
study  uses  a  1/3°  by  1/3°  grid.   Prior  simulations  of  the  1964  tsunami 
using  both  a  1/4°  by  1/6°  grid  and  a  1/5D  by  1/3°  grid  demonstrated 
negligible  differences  between  use  of  the  two  grid  spacings  (Houston  et 
al .  ,  1975a).   The  permanent  ground  deformation  used  as  input  is  also 
taken  from  Plafker  (1969).   The  numerical  grid  was  large  enough  to  cover 
a  substantial  portion  of  the  Pacific  Ocean  including  Wake  Island.   The 
numerical  model  calculations  at  Wake  Island  were  found  to  be  similar  to 
those  presented  by  Hwang  et  al.  (1972).   Figure  4  shows  a  comparison  of 
the  numerical  model  calculations  of  this  study  and  the  Wake  Island 
recording  of  the  1964  tsunami.   The  lack  of  agreement  in  the  trailing 
wave  region  may  be  due  to  local  effects  (e.g.  wave  trapping  by  the  Wake 
Island  seamount)  or  reflections  from  the  nearby  transmission  boundaries 
of  the  numerical  grid  (these  boundaries  are  not  perfectly  transmitting) . 

Figures  5  and  6  show  time-histories  of  the  1960  Chilean  and  1964 
Alaskan  tsunamis  in  deep  water  near  the  Hawaiian  Islands  calculated 
using  the  deep  ocean  propagation  model.   The  permanent  ground  defor- 
mation in  the  source  region  of  the  Chilean  earthquake  of  1960  is  taken 
from  Plafker  and  Savage  (1970).   Calculations  for  the  Alaskan  tsunami 
involve  a  grid  with  approximately  30,000  cells  and  the  calculations  for 
the  Chilean  tsunami  a  grid  with  approximately  90,000  cells. 

The  finite  element  numerical  model  is  a  time-harmonic  solution  of 
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the  boundary  value  problem.   The  response  of  a  group  of  islands  to  an 

arbitrary  tsunami  can  easily  be  determined  within  the  framework  of  a 

linearized  theory  using  the  theory  of  superposition.   For  example,  an 

arbitrary  tsunami  in  the  deep  ocean  can  be  Fourier  decomposed  as 

follows: 

b  (tO  =  f   b(w)  e"l(a)t+pfu)^  du  (2.36) 

where 

b  =  incident  wave  amplitude 
o  r 

b(w)  =  amplitude  of  frequency  component  tu 

p (u)  ~  phase  angle 

If  n(x,y,(j)   is  the  response  amplitude  at  any  point  (x,y)  along  the 

island  coasts  due  to  an  incident  plane  wave  of  unit  amplitude  and 

frequency  cj  ,  then  the  response  of  the  islands  to  the  arbitrary  tsunami 

time  history  b  ft)   is  given  by 
o 

f(x,y,t)  =  Re  [  /b(w)n(x,y,oO  e"1^^0^  *du]        (2.37) 

Therefore,  when  n(x,y,co)  is  known  for  all  u  ,  the  island  response  to 
an  arbitrary  tsunami  can  be  calculated.   Of  course,  it  is  not  feasible 
to  calculate  the  integrals  of  Equations  (2.56)  and  (2.37)  over  all 
frequencies.   Instead,  the  frequency  range  must  be  discretized  and  the 
integrals  replaced  by  sums  over  a  frequency  range  containing  most  of  the 
energy  of  the  tsunami. 

Equation  (2.36)  involves  a  Fourier  decomposition  of  a  time  series. 
This  decomposition  is  accomplished  for  the  time-histories  of  the  1960 
Chilean  and  1964  Alaskan  tsunamis  (Figures  5  and  6)  using  a  least 
squares  harmonic  fitting  procedure.   The  time-history  of  the  Alaskan 


tsunami  in  deep  water  was  decomposed  into  IS  components  with  periods 
ranging  from  14.5  minutes  to  the  time  length  of  the  record  (260 
minutes) .   The  variance  of  the  residual  (difference  between  the  actual 
record  and  a  recomposition  of  the  18  components)  was  approximately 
0.2  percent  of  the  variance  of  the  record.   Therefore,  virtually  all  the 
energy  of  the  wave  record  was  contained  in  the  18  components.   The 
time-history  of  the  Chilean  tsunami  in  deep  water  was  decomposed  into 
11  components  with  periods  ranging  from  15.5  minutes  to  the  time  length 
of  the  record  (170  minutes).   The  variance  of  the  residual  was  less  than 
0.1  percent  of  the  variance  of  the  record.   For  both  cases,  the  original 
time-history  and  a  time-history  constructed  from  a  recomposition  of  the 
components  were  virtual ly  indistinguishable. 

Equations  (2.56)  and  (2.57)  take  the  following  form  when  discre- 
tized : 

bo(t)  =  I  b(u)n)  e-i{V+p(V}  (2-5g) 

C(x,y,t)    =  Re[   E  b(u  )      n(x,y,ai  )    e"i{unt+pCwn5}]  (2.59) 

n=l  n 

where 

m  =  number  of  components, 

-  -   th 

to  =  frequency  or   n    component 

The  M"  )   term  is  determined  by  the  least- squares  harmonic 

fitting  procedure  and   n(x,y, u  )   is  determined  by  the  finite  element 

numerical  model  for  each  frequency  to   and  at  each  location  (x,y). 

Therefore,  the  time-history   5(x,y,t),  which  represents  the  response  of 

location  (x,y)  to  the  deepwater  tsunami  time-history  given  by  b  (t) , 


,s5 


can  be  calculated  for  any  location  along  the  coastline  of  the  Hawaiian 
Islands. 

Figures  7,  S,  and  9  present  comparisons  between  tide  gage  re- 
cordings and  numerical  model  calculations  of  the  1964  Alaskan  tsunami  at 
Kahului,  Maui;  Honolulu,  Oahu;  and  Hilo,  Hawaii,  respectively.   The 
largest  waves  recorded  at  each  of  these  sites  are  shown.   Waves 
arriving  at  later  times  are  all  much  smaller  than  those  shown.  When- 
ever tide  gage  limits  were  encountered,  the  recordings  were  linearly 
extended  (the  tide  gage  locations  are  shown  in  Figure  10) . 

The  wave  records  shown  in  these  figures  are  in  remarkable  agree- 
ment, especially  considering  the  fact  that  the  ground  displacement  of 
the  1964  earthquake  was  not  precisely  known.   Since  the  Hilo  breakwater 
was  not  included  in  the  numerical  model,  the  numerical  model  calcula- 
tions for  Hilo  are  probably  too  large;  however,  this  breakwater  was 
undoubtedly  highly  permeable  to  the  1964  tsunami.   The  numerical  results 
appear  to  be  too  large  in  Hilo  by  some  constant  factor,  since  the  tide 
gage  recording  and  the  numerical  model  calculations  have  the  same  form 
with  the  first  wave  crest  and  trough  having  approximately  the  same 
amplitude  and  being  proportionately  smaller  than  the  second  crest. 

Figure  11  shows  a  comparison  between  the  tide  gage  recording  of  the 
1960  Chilean  tsunami  at  Honolulu  and  the  numerical  model  calculations. 
Again  the  largest  waves  recorded  are  shown  (the  Honolulu  gage  was  the 
only  tide  gage  in  the  Hawaiian  Islands  not  destroyed  by  the  1960 
tsunami) .   The  permanent  ground  motion  of  the  1960  earthquake  is  not 
known  nearly  as  well  as  the  ground  motion  for  the  1964  earthquake. 
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Figurc  7.   1964  Tsunami  from  Alaska  Recorded  at  Kahului,  Maui 
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:igure  8.   1964  Tsunami  from  Alaska  Recorded  at  Honolulu,  Oahu. 
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Figure  9.   1961  Tsunami  from  Alaska  Recorded  at  fiiio,  Hawaii. 


Figure   10.      Gage    Location: 
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Hence,  the  agreement  indicated  in  Figure  11  is  quite  good  considering 
the  uncertainty  of  the  uplift  in  the  source  region. 

Tide  gages,  of  course,  do  not  record  tsunamis  perfectly  since  they 
are  nonlinear  devices  with  responses  which  depend  upon  both  the  period 
and  amplitude  of  the  disturbances  they  measure.   However,  simple 
calculations  based  upon  the  paper  of  Noye  (1970)  show  that  the  distortion 
is  small  for  the  tsunamis  shown  in  Figures  7  through  9  and  Figure  11. 

Figure  12  shows  the  response  amplitude  as  a  function  of  wave  period 
at  the  locations  of  tide  gages  in  Kahului,  Honolulu,  and  Hilo.   The 
response  versus  wave  period  is  greatly  dependent  upon  the  precise 
location  of  recording  for  bays  such  as  Hilo.   For  example,  the  Hilo 
response  (Figure  12)  at  the  tide  gage  location  shown  in  Figure  10  has 
a  single  large  resonant  peak  at  2S.9  minutes.   This  peak  occurs  at  the 
fundamental  period  of  oscillation  of  Hilo  Bay  and,  consequently,  all 
locations  in  the  bay  have  a  resonant  peak  at  this  incident  wave  period 
(the  amplitude  of  the  peak  depends  on  location  in  the  bay).   Other 
locations  in  Hilo  Bay  also  have  significant  peaks  at  lower  periods. 
Loomis  (1970)  calculated  a  free  oscillation  fundamental  period  (node 
assumed  at  bay  mouth)  of  26.2  minutes  for  Hilo  Bay  using  a  finite 
difference  grid  spacing  of  approximately  one-half  a  kilometer. 

Bernard  and  Vastano  (1977)  show  a  contour  map  of  a  normalized 
energy  parameter  versus  wave  period  and  location  around  the  island  of 
Hawaii.   The  map  shows  Hilo  displaying  an  energy  peak  from  10  to  16 
minutes  centered  at  approximately  12.5  minutes;  no  significant  peak  at 
28.9  minutes  is  observed.   The  difference  between  Bernard  and  Vastano 's 
calculations  and  the  results  presented  in  this  paper  is  attributed  to 
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the  poor  resolution  of  the  spatial  features  of  Hilo  Pay  by  their  finite 
difference  grid  spacing  of  5.5  kilometers.   Hi  Id  Bay  has  the  shape  of  a 
right  triangle  with  sides  that  are  approximately  9  kilometers  and  13.5 
kilometers  long.   The  grid  used  by  Bernard  and  Vastano  represented  Hilo 
Bay  as  a  right  triangle  with  sides  that  were  5.5  kilometers  and  11 
kilometers  long;  thus,  the  surface  area  of  Hilo  Bay  was  approximately 
one  half  its  actual  area.   This  distortion  of  the  sir.e  of  Hilo  Bay  and 
also,  perhaps,  the  distortion  of  the  bathymetry  of  the  bay  by  the  coarse 
grid  accounts  for  the  lower  response  periods  for  Hi lo  Bay  presented  by 
Bernard  and  Vastano. 

The  importance  of  accurately  representing  island  shapes  and  off- 
shore bathymetry  and  resolving  tsunami  wavelengths  can  be  illustrated  by 
considering  a  study  by  Adams  (1975)  that  used  the  numerical  model 
developed  by  Bernard  and  Vastano  (19""")  to  interpolate  between  his- 
torical measurements  of  tsunami  elevations  produced  by  the  1946  Aleutian 
tsunami.   A  single  Gaussian  wave  crest  was  used  as  input  for  the 
numerical  model.   If  a  resolution  of  four  grid  cells  per  crest  wave- 
length is  maintained,  this  numerical  model  cannot  propagate  waves  into 
water  with  a  depth  less  than  about  500  feet.   The  main  processes  in- 
volved in  the  transformation  of  the  form  of  a  tsunami  probably  occur  in 
the  region  over  the  continental  shelf  (depths  less  than  600  feet). 
Consequently,  the  elevations  calculated  by  Adams  were  small  and  did  not 
display  the  great  variation  along  coastlines  known  to  occur  during 
historical  tsunamis  in  the  Hawaiian  Islands. 

Adams  included  the  effects  of  propagation  over  continental  shelves 
by  multiplying  the  numerical  model  elevations  by  a  factor  (different 
for  each  island)  that  forced  the  best  agreement  with  the  historical 
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data.   However,  even  with  artiface,  agreement  with  historical  data  was 
poor.   For  example,  on  the  northeast  coast  of  Oahu,  the  numerical  model 
elevations  multiplied  by  the  factor  resulted  in  predicted  elevations 
ranging  from  14  to  21.4  feet.   Loomis  (1976)  presents  36  historical 
elevations  for  the  1946  Aleutian  tsunami  on  the  northeast  coast  of  Oahu 
and  only  7  even  fall  within  this  range  of  elevations.   The  historical 
elevations  actually  varied  from  1  to  37  feet  on  this  coast. 


CHAPTER  III 

TSUNAMI  INTERACTION  WITH  CONTINENTAL  COASTLINES 

Introduction 

Experimental  investigations  of  long  wave  interactions  with 
continental  coastlines  have  generally  been  one-dimensional  and  have 
involved  idealized  conditions.   Flume  tests  by  Savage  [1358)  and 
Saville  (1956)  measured  runup  of  periodic  waves  on  constant  slope 
beaches  and  Mail  and  Watts  (1953),  Kaplan  (1935),  and  Kishi  and 
Saeki  (1966)  measured  runup  of  solitary  waves  :>n  constant  slope 
beaches.   All  of  these  flume  experiments  used  waves  much  steeper 
than  typical  tsunamis  due  to  the  difficulties  involved  in  the 
generation  of  small  amplitude  long  waves. 

Many  theoretical  studies  of  one-dimensional  propagation  of 
weakly  dispersive  long  waves  have  been  performed  in  recent,  years. 
Tsunamis  become  weakly  dispersive  (amplitude  and  frequency  dispersive) 
long  waves  at  some  point  during  their  propagation,  although  (as  was 
the  case  for  the  Hawaiian  Islands)  they  may  not  be  weakly  dispersive 
for  a  time  period  sufficiently  long  for  related  effects  to  develop. 
The  proper  equations  for  the  investigation  of  weakly  dispersive 
long  wave  propagation  over  a  steeply  sloping  bottom  were  first 
derived  by  Mei  and  LeMehaute  (1966)  and  in  a  slightly  different 
form  by  Peregrine  (1967)  who  also  gave  the  first  numerical  solution 
in  the  case  of  a  solitary  wave  using  an  implicit  finite  difference 
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scheme.   A  more  complete  solution  was  later  given  by  Madsen  and  Mei 
(1969)  who  transformed  the  equations  to  characteristic  form  before 
they  were  solved  numerically.  Svendsen  (1974)  studied  weakly  dis- 
persive wave  propagation  over  a  gently  sloping  bottom.   Hammack 
(1972)  made  comparisons  of  long  wave  generation  by  bottom  displace- 
ments in  the  laboratory  and  numerical  solutions  of  the  Korteweg- 
DeVries  equation  using  the  numerical  approach  of  Peregrine  (1966). 

Several  studies  of  one-dimensional  tsunami  propagation  over 
continental  shelves  of  constant  slope  have  been  performed.   Keller 
and  Keller  (1964)  solved  nonlinear  long  wave  equations  numerically 
for  a  constant  slope  bathymetry.   Takahasi  (1964)  studied  tsunami 
propagation  over  a  simple  bathymetry  through  analytic  solutions  of 
the  linear  long  wave  equation.   Carrier  (1966)  assumed  a  tsunami 
was  generated  by  a  point  source  and  studied  nonlinear  propagation 
over  a  constant  slope  bathymetry.   Heitner  and  Housner  (1970)  used 
a  finite  element  numerical  method  to  solve  nonlinear  equations  for 
propagation  up  constant'  slope  beaches.   More  recently,  Mader  (1974) 
studied  tsunami  propagation  over  linear  slopes  using  a  finite 
difference  solution  of  the  complete  Navier-Stokes  equations  (not 
vertically  integrated). 

Occasional  two-dimensional  studies  have  been  performed  for 
tsunami  propagation  in  ncarshore  regions.   Grace  (1969)  studied 
tsunamis  affecting  the  reef  runway  of  the  Honolulu  International 
Airport  using  a  hydraulic  model.   However,  the  wavemakers  used  were 
only  a  small  fraction  of  a  tsunami  wavelength  from  the  region  to  be 
studied  and  waves  reflected  by  land  were  reflected  again  almost 
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immediately  by  the  wavemakers.  This  is  a  common  problem  in  hydraulic 
model  tests  involving  long  waves.   The  wavemakers  can  be  moved  a 
few  wavelengths  from  the  region  of  interest  by  reducing  the  model 
scale;  however,  scale  effects  may  then  become  very  significant. 
Wave  absorber  screens  may  also  be  placed  in  front  of  the  wavemakers; 
however,  it  is  very  difficult  to  absorb  very  long  waves  and  such  an 
approach  is  often  not  very  helpful.   Similar  problems  with  re- 
reflected  waves  can  occur  in  numerical  studies.   The  calculations 
of  Hwang  and  Divoky  (1975)  were  in  poor  agreement  with  the  tide  gage 
recording  of  the  1964  tsunami  at  Hilo  not  only  for  the  reasons 
presented  in  Chapter  II,  but  also  because  their  algorithm  allowed 
waves  to  leave  the  numerical  grid  efficiently  only  for  waves  approach- 
ing the  input  boundary  normally.   Since  Hilo  Bay  is  triangularly 
shaped,  the  reflected  waves  approached  the  input  boundary  at  various 
angles  and  were  partially  re-reflected  toward  Hilo. 

Governing  Equations 

The  governing  equations  for  long  wave  propagation  over  the 
continental  shelf  of  the  west  coast  of  the  United  States  are  Equations 
(2.5)  and  (2.4).   Thus,  the  long  wave  assumption  is  made  once  again 
with  vertical  velocities  and  accelerations  neglected  relative  to 
horizontal  components.   The  continuity  equation  is  given  by  the 
equation 

~   +  V  •  qCn  +  h)  =  0  (3.1) 

and  the  momentum  equation  is 


8q   ->    -v  2->- 

yf  +  q  •  Vq  =  pF  -  Vp  +  pV  q  (3.2) 
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If  the  only  external  force  is  gravity  and  the  pressure  is  hydro- 
static, then  liquation  (3.2)  becomes 

p  S- +  q  •  Vci  =  psVri  +  ^iV"ci  (3-^) 

c)t 

The  west  coast  of  the  United  States  has  a  much  broader  continental 
shelf  than  that  of  the  Hawaiian  Islands;  thus,  the  nonlinear  advective 
and  bottom  stress  effects  may  be  more  significant,  and  the  nonlinear 
and  viscous  terms  of  Equation  (3.3)  must  be  retained. 

The  dynamic  viscosity  u   in  the  viscous  term  of  Equation  (3.3) 
is  replaced  by  an  eddy  viscosity  constant   e  for  the  case  of 
turbulent  flow.   Using  this  approach  the  viscous  term  can  be  expressed 
as 

*,  c^4  ♦  ^4)  +  b2  M  C3.4) 

1   3x"    3y         ^z 

where  the  coefficients   e.  and  e_  represent  coefficients  of 

horizontal  and  vertical  eddy  viscosity,  respectively.   The  term 

2     2 

3x    3y 

82 
is  usuallv  small  in  comparison  with  the  e     — §■   term  (Dronkers, 

1964)  and  will  be  neglected  in  this  study.     If  e,  — |  is 

9z 
integrated  with  respect  to  z  from  the  bottom  to  the  surface 

and  divided  by  the  distance  from  bottom  to  surface  n  +  h  ,  we  find 


[(ItK-.k  -  (It)  .  -  J  V-V 


n+h    L  9z  z=-h    3z   z  =  n 
The  terms   £:_(-—)   _  .   and   E2('^z=n  &re  coraPonents  of  the 
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tangential  stresses  at  the  sea  bottom  and  surface,  respectively. 
The  shear  stress  at  the  sea  surface  is  quite  important  for  storm 
surge  problems  where  wind  velocities  are  very  large  but  can  be 
neglected  relative  to  the  bottom  stress  for  tsunamis. 

From  experiments  and  theoretical  considerations  for  one- 
dimensional  flow  in  a  river,  it  has  been  shown  (Dronkers,  1964),  that 
the  bottom  stress   t,   can  be  expressed  as  follows 

T  =l£SHM  (3.6) 

b    c2 

where  C   is  Chezy's  coefficient  and  is  given  by  the  expression 

1,486  hl/6  (37) 

n 
m 

where  n   is  Manning's  n  .   Thus,  Equation  (3.3)  may  be  written 
m 

-*-,  2   2,1/2 
3q   •*"   ^            n    pgq(u  +v  )  f~   R1 

p  —x   +  q  •  Vq  =  pgVr.  -  ~^ — ' ' U-6J 

CT(n+h) 

Equation  (3.8)  contains  the  nonlinear  advective  term  but  neglects 
frequency  dispersion  since  vertical  accelerations  have  been  neglected, 
This  equation  adequately  describes  wave  propagation  for  Ursell 
numbers  much  greater  than  unity  (Ursell,  1953).   The  Ursell  number 

is  proportional  to  the  ratio  of  nonlinear  effects  to  linear  effects 

nL2 
(i.e.,  frequency  dispersion)  and  is  equal  to  — =—  ,  where  L  is  a 

1/2         h 
typical  wavelength.  Since  L=(gh)   T  for  long  waves,  where  T  is 

T  2 
the  wave  period,  an  Ursell  number  may  be  chosen  as  gn(^) 

The  numerical  model  described  in  the  next  section  uses  Equation 

(3.8)  to  propagate  waves  from  depths  of  500  meters  to  shore.  These 
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waves  have  periods  of  the  order  of  1800  seconds  and  initial  wave 

T  2 
amplitudes  of  the  order  of  a  foot.   Thus,  the  Ursell  number,   gnGr) 

is  aproximately  50  and  increases  as  the  wave  enters  shallow  water 

with  h  decreasing  and   n  increasing. 

Mader  (1974)  found  that  the  one-dimensional  form  of  Equation 

(3.8),  without  the  viscous  term,  produced  very  similar  results  to  a 

one-dimensional  N'avier-Stokes  equation  (viscosity  neglected)  for 

propagation  of  long  wavelength  tsunamis  over  linear  continental 

slopes.   The  pressure  term  in  Mader' s  solution  was  not  hydrostatic; 

thus,  frequency  dispersion  was  permitted. 

Finite  Difference  Model 


To  solve  Equations  (5.1)  and  (5.8)  numerically,  the  differential 
equations  are  ciscretized  and  replaced  by  a  system  of  finite  difference 
equations  using  central  differences  on  a  space-staggered  grid.   The 
space-staggered  scheme  describes  velocities,  water  levels,  and  depths 
at  different  grid  points.   This  scheme,  first  used  by  Platzman  (195S) 
for  lake  surge  problems,  has  the  advantage  that  in  the  equation  for  the 
variable  operated  upon  in  time,  there  is  a  centrally  located  spatial 
derivative  for  the  linear  term.   Leendertse's  (1967)  implicit-explicit 
multi-operational  method  is  employed  in  determining  the  solution 
for  q  and   n  as  functions  of  time.   The  stability  and  convergence 
of  this  method  is  described  in  detail  by  Leendert.se  (1967)  .   This 
method  for  solving  finite-difference  equations  has  been  applied  in 
studies  analyzing  the  tidal  hydraulics  of  harbors  and  inlets  (Raney, 
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1976).   Houston  and  Garcia  (1975)  used  the  method  to  calculate  the 
decay  of  a  tsunami  as  it  enters  and  spreads  throughout  San  Francisco 
Bay. 

The  nearshore  finite-difference  numerical  model  uses  the  time- 
history  calculated  by  the  generation  and  deep  ocean  propagation 
model  described  in  Chapter  II  as  input  in  the  following  manner.   A 
tsunami  is  generated  in  the  Aleutian-Alaskan  area  or  the  west  coast 
of  South  America  and  propagated  across  the  deep  ocean  to  a  500- 
meter  depth  off  the  west  coast  of  the  United  States.   Waveforms 
calculated  at  this  depth  by  the  deep  ocean  numerical  model  are 
recorded  all  along  the  west  coast.   These  waveforms  are  then  used 
as  input  to  the  nearshore  numerical  model,  which  propagates  the 
tsunamis  from  the  500-meter  depths  across  the  continental  slope  and 
shelf  to  shore. 

Figure  13  shows  outlines  of  the  four  numerical  grids  used  to 
cover  the  west  coast.   The  grids  have  square  grid  cells  two  miles 
on  a  side  (southern  California  and  some  of  San  Francisco  Bay  were 
not  covered  by  numerical  grids  because  elevation  predictions  were 
made  in  previous  studies  (Houston  and  Garcia,  1974;  and  Houston  and 
Garcia,  1975)  using  simple  one-dimensional  analytical  solutions  for 
propagation  over  the  continental  shelf).   The  offshore  bathymetry 
was  modeled  from  the  500-meter  contour  to  shore.   Beyond  the  500- 
meter  contour,  the  ocean  was  assumed  to  have  a  constant  500-meter 
depth.   The  input  boundary  of  each  grid  was  located  approximately 
one  and  one-half  wavelengths  of  a  30-minutc  wave  from  the  shore. 
Therefore,  at  least  three  typical  waves  could  arrive  at  the  shore 
before  waves  reflected  from  the  input  boundary  became  a  problem. 
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Figure  15.   Finite  Difference  Grid  Locations, 
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The  input  boundaries  of  the  grids  were  oriented  approximately 
parallel  to  the  shoreline  since  refraction  will  bend  the  wavefronts 
to  such  an  orientation  before  they  reach  the  500-meter  contour 
(Wilson  and  Torum,  1968).   Lateral  boundaries  of  the  grids  were 
taken  to  be  impermeable  vertical  walls. 

Model  Verification 

The  finite  difference  nearshcre  numerical  model  was  also  veri- 
fied by  numerical  simulations  of  the  1964  Alaskan  tsunami.   The 
deep  ocean  propagation  model  discussed  in  Chapter  II  was  used  to 
generate  the  1964  tsunami  and  propagate  it  across  the  deep  ocean. 
Figure  14  shows  surface  elevation  contours  3-1/2  hours  after  the 
1964  Alaskan  earthquake,  as  calculated  by  the  deep-ocean  propagation 
numerical  model.   This  figure  illustrates  the  concentration  of  energy 
on  the  northern  California,  Oregon,  and  Washington  coasts  due  to  the 
directional  radiation  of  energy  from  the  source  region.   The  upper 
section  of  Figure  15  shows  a  time  history  of  the  1964  Alaskan  tsunami 
calculated  by  the  deep-ocean  model  near  Crescent  City,  California,  in 
a  water  depth  of  500  r.eters.   This  waveform  was  used  as  input  to  the 
finite  difference  nearshore  model. 

The  lower  section  of  Figure  15  shows  a  comparison  between  a  tide 
gage  recording  of  the  1964  tsunami  at  Crescent  City  and  the  nearshore 
model  simulation  of  this  tsunami.   The  numerical  model  calculations 
agree  well  with  the  tide  gage  recording;  periods,  phases,  and 
amplitudes  are  accurately  reproduced.   The  main  disagreement  occurs 
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Figure  15.   1964  Tsunami  from  Alaska  Recorded  at  Crescent  City,  California 
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for  the  amplitude  of  the  third  wave.   With  regard  to  this  disagree- 
ment, it  should  be  noted  that  there  is  also  a  lack  of  agreement 
on  the  amplitudes  of  the  actual  historical  third  and  fourth  wave 
crests.   Figure  15  shows  a  reconsti action  of  the  1964  tsunami  at 
Crescent  City  inferred  by  Wilson  and  Torum  (1963)  from  the  prototype 
gage  record  and  later  survey  measurements  by  Magoon  (1965).   The 
elevation  of  the  fourth  wave  (not  shown)  was  estimated  by  Wilson 
and  Torum  to  have  been  a  little  less  than  14  feet  above  mean  lower 
low  water  (MLLW)  (the  zero  elevation  shown  in  Figure  15  is  MLLW 
datum);  however,  Wiegel  (1965)  estimated  that  the  elevation  of  the 
third  wave  was  approximately  16  feet  above  MLLW  and  that  the  fourth 
wave  attained  the  highest  elevation  of  18  or  19  feet  above  MLLW  at  the 
tide  gage  location.   The  numerical  mode]  predicted  a  fourth  wave 
elevation  approximately  the  same  as  that  estimated  by  Wilson  and 
Torum  (14  feet  above  MLLW) .  However,  reflections  from  the  input 
boundary  were  probably  growing  in  importance  in  the  numerical  model 
calculations  during  the  arrival  of  the  fourth  wave. 

There  may  be  several  reasons  why  the  maximum  wave  elevation 
predicted  by  the  numerical  model  does  not  agree  with  the  historical 
record  as  well  as  the  elevations  of  the  initial  waves.   First,  the 
later  waves  in  the  deepwater  waveform  are  probably  waves  which  have 
been  reflected  from  land  areas  in  the  source  region.   The  fairly 
large  grid  cell  spacing  used  in  the  deep  ocean  propagation  model 
to  represent  the  land-water  boundary  in  this  region  probably  distorts 
the  reflected  waveform.   Total  reflection  of  waves  in  this  region 
also  is  not  completely  realistic.   Furthermore,  the  nearshore  grid  may 
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not  have  fine  enough  cells  to  completely  model  the  resonance 
phenomenon  leading  to  the  large  third  or  fourth  wave. 

Since  1964,  there  has  been  considerable  speculation  concerning 
the  reasons  that  the  effects  of  the  1964  tsunami  were  so  great  at 
Crescent  City.  .The  finite  difference  numerical  models  show  that 
both  the  directional  radiation  of  energy  from  the  source  region 
(see  Figure  14)  and  a  local  resonance  caused  the  relatively  large 
elevations  at  Crescent  City.   Actually,  large  elevations  at  Crescent 
City  were  not  unique  since  directional  radiation  of  energy  from  the 
source  also  caused  large  elevations  along  the  Oregon  and  northern 
California  coasts.   Runup  10  to  15  feet  above  the  high  tide  level 
occurred  all  along  the  Oregon  coast  south  of  the  Columbia  River 
(runup  at  Crescent  City  was  approximately  15  feet  above  the  high 
tide  level).   However,  the  Oregon  coast  is  very  sparsely  populated 
and  there  were  few  damage  reports.   The  severity  of  structural 
damage  at  Crescent  City,  which  has  a  large  logging  industry,  also 
was  apparently  due  to  the  impact  of  logs  carried  by  the  tsunami 
(Wilson  and  Torum,  1968) .   The  nearshore  numerical  model  indicates 
that  the  resonant  effects  at  Crescent  City  were  fairly  local, 
extending  over  2  to  4  miles  of  coastline.   This  behavior  contradicts 
Wilson  and  Torum' s  (1968)  speculation  that  a  bowl-shaped  section  of 
the  continental  shelf  with  a  diameter  of  approximately  50  miles 
experienced  a  resonant  oscillation.   Historical  data  support  the 
numerical  model  calculations.   At  the  mouth  of  the  Klamath  River, 
approximately  15  miles  south  of  Crescent  City,  elevations  observed 
during  the  1964  tsunami  were  only  2  to  3  feet  above  normal  high 
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tide.  Wilson  and  Torum's  hypothesis  of  a  shelf  oscillation  would 
predict  elevations  at  the  mouth  of  the  Klamath  River  greater  than 
those  that  occurred  at  Crescent  City. 

Figure  16  shows  a  comparison  between  a  tide  gage  recording  of 
the  1964  tsunami  at  Avila  Beach,  California,  and  the  nearshore 
numerical  model  calculations.   The  elevations  recorded  by  the  Avila 
Beach  tide  gage  were  larger  than  those  recorded  at  any  tide  gage  on 
the  west  coast  except  the  Crescent  City  gage.   The  historical  and 
calculated  waveforms  shown  in  Figure  16  are  in  good  general  agreement. 
The  tide  gage  record  obviously  has  higher  frequency  components  which 
are  not  predicted  by  the  numerical  model;  these  components  may  be  local 
oscillations  of  water  areas  which  were  too  small  to  be  accurately 
represented  by  the  numerical  grid.   Important  features  such  as  the  wave 
amplitudes  are  in  good  agreement. 

The  waveforms  recorded  for  the  1964  tsunami  by  the  Crescent 
City  and  Avila  Beach  tide  gages  have  elevations  larger  than  those 
recorded  by  any  other  gage  on  the  west  coast.   Other  tide  gage 
recordings  of  this  tsunami  on  the  west  coast  (except  at  Astoria, 
Oregon)  were  in  areas  not  covered  by  the  grids  shown  in  Figure  13. 
A  comparison  with  the  tide  gage  at  Astoria,  Oregon,  was  not  made 
because  Astoria  is  approximately  12  miles  away  from  the  coast  in 
the  estuary  of  the  Columbia  River,  and  amplitudes  of  the  1964 
tsunami  were  small  at  Astoria. 
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CHAPTER  IV 

ELEVATION  PREDICTIONS 

Hawaiian  Island  Predictions 

Introduction 

The  Hawaiian  Islands,  a  chain  of  eight  islands  as  shown  in 
Figure  17,  have  a  history  of  destructive  tsunamis  generated  both  in 
distant  areas  and  locally.   The  earliest  recording  of  a  severe  tsunami 
in  the  Hawaiian  Islands  was  in  1S37  when  a  tsunami  from  Chile  reached 
an  elevation  of  20  feet  at  Hilo,  Hawaii,  and  killed  46  people  in  the 
Kau  District  of  the  big  island  of  Hawaii.   Prior  to  1S37,  a  number 
of  severe  tsunamis  undoubtedly  reached  the  islands  but,  unfortu- 
nately no  detailed  records  were  kept.   Since  1S37,  there  have  been  16 
tsunamis  that  have  caused  significant  damage. 

Most  of  the  destructive  tsunamis  in  the  Hawaiian  Islands  have 
been  generated  along  the  coast  of  South  America,  the  Aleutian  Islands, 
and  the  Kamchatkan  Peninsula  of  the  Soviet  Union.   Approximately  one- 
fourth  of  all  the  tsunamis  recorded  in  the  Hawaiian  Islands  have 
originated  along  the  coast  of  South  America,  while  more  than  one- 
half  have  originated  in  the  Kuril-Kamchatka-Aleutian  region  of  the 
north  and  northwestern  Pacific.   Tsunamis  generated  by  local  seismic 
events  have  caused  large  runup  in  the  islands,  especially  on  the 
southeast  coast  of  the  big  island  of  Hawaii.   The  1868  tsunami 
produced  the  largest  waves  of  record  in  the  Hawaiian  Islands  with 
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60-foot  waves  reported  on  the  South  Puna  coast  of  the  island  of 
Hawaii.   The  most  recent  major  tsunami  in  Hawaiian  history  occurred 
on  November  29,  1975,  when  waves  generated  by  an  earthquake  with 
an  epicenter  on  the  South  Puna  Coast  may  have  reached  elevations  as 
great  as  45  feet  along  the  southeast  coast.   The  most  destructive- 
tsunami  to  ever  hit  the  islands  in  terms  of  both  loss  of  life  and 
property  destruction  was  the  Great  Aleutian  tsunami  of  1946,  which 
killed  175  people  and  produced  waves  over  55  feet  in  elevation.   Hilo 
incurred  $26  million  in  property  damage  attributable  to  this  tsunami. 
The  1960  Chilean  tsunami  is  the  most  recent  distantly  generated 
tsunami  that  produced  major  effects  in  the  Hawaiian  Islands.   Sixty- 
one  lives,  all  at  Hilo,  were  claimed  by  the  tsunami.   Damage  through- 
out the  state  was  estimated  to  be  $25.5  million.   Inspection  of  the 
damage  at  Hilo  revealed  much  evidence  of  the  tremendous  forces 
developed  by  the  waves.   Twenty-ton  boulders  had  beer,  moved  hundreds 
of  feet,  asphaltic  concrete  pavements  were  peeled  from  their  sub- 
base,  and  hundreds  of  automobiles  were  moved  and  crushed. 

The  following  sections  discuss  the  methods  used  to  obtain 
sufficient  data  to  predict  runup  elevations  for  the  entire  coastline 
of  the  Hawaiian  Islands.   The  finite  element  numerical  model  discussed 
in  Chapter  II  is  used  to  interpolate  between  local  historical  data 
for  the  years  of  accurate  survey  measurements  from  1946  through  1977. 
Such  a  method  of  interpolation  is  needed  so  that  elevations  can  be 
predicted  for  the  entire  coastline  even  for  locations  having  no 
actual  historical  data.   A  method  is  discussed  herein  which  allows 
reconstruction  of  historical  data  for  the  period  of  time  from  1857  to 
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1946  when  elevations  for  the  period  from  1946  through  1977  are  known. 
Such  an  approach  is  necessary  in  order  to  have  sufficient  data  for 
accurate  elevation  predictions  since  it  is  shown  herein  that  reason- 
able predictions  cannot  be  made  solely  based  upon  known  historical 
data  for  the  short  time  span  from  1946  through  1977. 
Interpolation  of  Recent  Historical  Pat a 

Historical  data  of  tsunami  activity  in  the  Hawaiian  Islands 
are,  of  course,  often  limited  to  certain  locations.   Information  on 
tsunami  activity  in  the  islands  prior  to  the  1946  tsunami  is  concentrated 
in  Hilo  (Hawaii)  and  to  a  lesser  extent  in  Honolulu  (Oahu)  and 
Kahului  and  Lahaina  (Maui).   Even  data  for  tsunamis  occurring  during 
the  period  of  accurate  survey  measurements  from  194  6  through  the 
present  are  absent  or  fragmentary  (i.e.,  data  exist  only  for  certain 
of  the  events)  for  most  of  the  coastline  of  the  islands.   Therefore, 
it  is  necessary  to  rely  on  more  than  just  available  historical  data 
to  determine  tsunami  occurrence  frequencies. 

The  numerical  model  is  used  to  fill  in  historical  data  gaps  for 
tsunamis  from  1946  through  1977  by  providing  relative  responses  of 
the  Hawaiian  Islands  to  tsunamis.   For  example,  although  the  deepwater 
waveform  of  a  tsunami  such  as  that  of  1946  is  not  known,  the  direction 
of  approach  of  this  tsunami  and  its  range  of  wave  periods  are  known. 
The  average  response  of  the  islands  to  a  tsunami  similar  to  the  1946 
tsunami  is  determined  by  inputting  sinusoidal  waves  of  unit  amplitude 
from  the  direction  of  the  1946  tsunami  into  the  numerical  model  over 
a  band  of  wave  periods.   If  historical  data  exist  at  one  location, 
wave  elevations  at  a  nearby  location  for  which  historical  data  do 
not  exist  are  determined  by  multiplying  the  historical  data  at  the 
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first  location  by  the  ratio  of  the  response  calculated  by  the  numer- 
ical model  at  the  second  location  and  the  response  calculated  at  the 
first.   The  numerical  model  takes  into  account  the  major  processes 
that  would  cause  different  relative  wave  elevations  at  the  two 
locations.  That  is,  the  model  calculates  shoaling,  refraction, 
diffraction,  reflection,  resonance,  shielding  of  the  backside  of  an 
island  by  the  front  side,  and  interactions  between  islands.   The 
historical  data  account  for  the  unknown  absolute  heights  of  the 
deepwater  tsunamis.   The  "correction  factors"  provided  by  the 
historical  data  are  local  instead  of  global  (single  factor  for  all 
islands)  because  there  are  certain  local  factors  that  cannot  be 
determined  by  the  numerical  model.   For  example,  coral  reefs  extend  over 
parts  of  some  of  the  islands  and  are  known,  to  protect  against  tsunami 
attack  (e.g.,  the  extensive  reefs  in  Kaneohe  Bay,  Oahu,  which  are 
powerful  tsunami  dissipators) .   The  effects  of  the  reefs  on  tsunamis 
are  implicitly  contained  in  the  historical  data  itself.   The  historical 
data  from  1946  through  the  present  were  taken  from  the  most  recent 
and  complete  compilation  of  these  data  (Loomis,  1976)  . 

As  an  example  of  the  use  of  the  numerical  model  to  fill  in 
historical  data  gaps  for  tsunamis  from  194  6  through  the  present, 
consider  the  1946  tsunami  recorded  on  the  island  of  Lanai  only  at 
Kaumalapua  Harbor  and  Manele  Bay.   Wave  elevations  can  be  calculated 
for  the  1946  tsunami  at  other  locations  on  the  coastline  of  Lanai  by 
multiplying  these  historical  values  by  the  ratios  of  the  responses 
calculated  by  the  numerical  model  at  locations  of   interest  and  the 
responses  calculated  at  Kaumalapua  Harbor  and  Manele  Bay.   In  this 


way,  elevations  for  the  1946  tsunami  can  be  determined  for  all  of 

Lanai . 

Time  Period  Analysis 

It  is  necessary  to  use  historical  data  of  tsunami  occurrence 
covering  the  greatest  possible  time  span  to  properly  determine 
frequency  of  occurrence  of  tsunami  elevations  for  the  Hawaiian 
Islands.   Tsunami  activity  has  not  been  uniform  in  these  islands 
during  recorded  history.  For  example,  the  two  largest  and  four  of 
the  ten  largest  tsunamis  striking  Hilo  from  1857  through  1977  occurred 
during  the  15-year  period  from  1946  through  1960.   Two  of  the  tsunamis 
from  1946  through  1960  originated  in  the  Aleutian  Islands,  one  in 
Kamchatka  and  one  in  Chile.  However,  six  of  the  ten  largest  tsunamis 
occurred  during  the  109-year  period  from  1857  to  1945  with  three 
originating  in  Chile,  two  in  Kamchatka  and  one  in  Hawaii.   Therefore, 
both  the  frequency  of  occurrence  and  place  of  origin  of  tsunamis 
have  been  remarkably  variable.   Any  study  (e.g.,  Towill  Corporation, 
1975)  basing  frequency  calculations  on  a  short  time  span  that  includes 
the  period  from  1946  through  1960  will  predict  a  significantly  more 
frequent  occurrence  of  large  tsunamis  than  is  warranted  by  historical 
data  from  1857  through  1977. 

The  errors  introduced  in  frequency  of  occurrence  calculations 
by  consideration  only  of  a  short  time  period  that  includes  the 
unrepresentative  years  from  1946  through  1960  will  be  greater  than 
the  errors  resulting  from  possible  observational  inaccuracies  of  the 
19th  Century.   This  is  easily  demonstrated  for  Hilo,  Hawaii,  which 
is  the  location  in  the  Hawaiian  Islands  having  the  most  complete 
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data  of  tsunami  activity.   Using  the  data  compiled  by  Cox  (1964)  for 
the  ten  largest  tsunamis  in  Hilo  from  1857  through  197"  (his  data 
were  actually  through  1964,  but  no  large  tsunamis  have  occurred  in 
Hilo  since  1964),  a  least-squares  fit  of  the  data  was  made  employing 
the  logarithmic  frequency  distribution  discussed  later  in  this 
chapter.   The  100-year  tsunami  elevation  for  Hilo  calculated  using 
this  method  was  27.3  feet.  If  the  elevations  recorded  during  the 
nineteenth  centry  (five  of  the  ten  largest  tsunamis),  are  increased 
by  50  percent  (which  is  much  larger  than  possible  observational 
error),  the  100-year  elevation  is  found  to  be  50.4  feet.   Similar 
calculations  just  based  upon  data  taken  during  the  period  of  accurate 
survey  measurements  in  Hilo  from  194  6  through  1977  yields  a  100-year 
elevation  of  44.2  feet.   Since  the  largest  elevation  in  Cox's  data 
for  the  141-year  period  from  1S57  through  1977  was  just  28  feet,  the 
100-year  elevation  of  44.2  feet  is  obviously  too  large.   Similar 
overestimates  will  occur  at  any  location  in  the  Hawaiian  Islands  if 
an  analysis  is  based  upon  a  short  period  of  time  that  includes  the 
years  1946  through  1964  since  the  exceptionally  frequent  occurrence 
of  major  tsunamis  in  Hilo  from  1946  through  1964  is  a  property  of 
the  unusual  activity  of  tsunami  generation  areas  and  not  of  special 
properties  of  Hilo. 

The  unusual  tsunamigenic  activity  in  generation  regions  in 
recent  years  is  reflected  in  the  most  up-to-date  catalog  of  tsunami 
occurrence  in  the  circumpacific  area  complied  by  Soloviev  and  Go 
(1969).   This  catalog  shows,  for  example,  that  the  three  greatest 
intensity  tsunamis  generated  in  the  Aleutian-Alaskan  region  since 
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1788  occurred  in  1946,  1957,  and  1964.   The  catalog  also  lists  the 
1960  Chilean  tsunami  as  the  greatest  intensity  tsunami  generated  in 
South  America  in  recorded  history  [416-year  period  from  1562  through 
1977)  .   Reporting  of  tsunamis  on  the  west  coast  of  South  America  is 
quite  good  with, 105  tsunamis  reported  during  the  416-year  period. 
Kanamori  (1977]  also  notes  the  unusual  seismic  activity  that  occurred 
during  a  relatively  short  time  span  in  recent  years.   He  shows  that 
the  four  earthquakes  having  the  greatest  seismic  moments  (true 
measure  of  magnitude  for  great  earthquakes)  in  this  century  occurred 
during  a  thirteen  year  period  from  1952  through  1964  (1960  Chilean, 
1964  Alaskan,  1957  Aleutian,  and  1952  Kamchatkan  earthquakes). 

Evidence  that  tsunamis  were  unusually  active  in  the  Hawaiian 
Islands  from  1946  through  1964  is  also  apparent  from  historical  data 
at  Kahului  and  Lahaina,  Maui.   Kahului,  Lahaina,  iiilo,  and  Honolulu 
are  the  only  locations  in  the  islands  with  historical  data  extending 
to  1837.   Four  of  the  five  largest  tsunami  elevations  recorded  at 
Kahului  occurred  during  the  years  from  1946  through  1964.   Three  of 
the  four  largest  tsunami  elevations  recorded  at  Lahaina  occurred 
during  the  same  period.   The  historical  data  for  Honolulu  is  not 
very  useful  because  the  commonly  reported  elevations  for  the  1837 
and  1841  tsunamis  were,  in  fact,  drops  in  water  level  (Pararas- 
Carayannis,  1977)  and  there  are  no  reports  for  the  1868  and  1877 
tsunamis  from  Chile. 
Reconstruction  of  Historical  Data 

In  order  to  have  data  covering  a  sufficiently  long  time  span  to 
make  reasonable  predictions  of  tsunami  elevations  for  the  Hawaiian 
Islands,  it  is  necessary  to  reconstruct  historical  data  for  the 
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years  prior  to  1946.   To  do  this,  it  is  assumed  in  this  study  that 
tsunamis  generated  in  a  single  source  region  (Kamchatka  or  Chile  but 
not  the  Aleutian-Alaskan  region)  approach  the  islands  from  approximately 
the  same  direction  and  have  energy  lying  within  the  same  band  of 
wave  periods.   The  difference  in  wave  elevations  at  the  shoreline  in 
the  Hawaiian  Islands  produced  by  tsunamis  generated  at  different 
times  in  the  same  region  is  attributed  mainly  to  differences  in 
deepwater  wave  amplitudes.   For  example,  the  1841  tsunami  from 
Kamchatka  produced  a  wave  elevation  in  Hilo  that  was  approximately 
25  percent  greater  than  that  of  the  1952  tsunami  from  Kamchatka.   It 
is  assumed  that  these  two  tsunamis  had  the  same  relative  magnitudes 
throughout  the  islands.  Since  the  finite  element  model  allows  a 
determination  of  elevations  for  the  1952  tsunami  along  the  entire 
coastline  of  the  islands,  the  elevations  of  the  1S41  tsunami  and  all 
ether  Kamchatkan  tsunamis  occurring  during  the  period  of  record  in 
Hilo  from  1857  to  1946  can  be  reconstructed  for  the  entire  coastline. 
Similarly,  Chilean  tsunami  elevations  can  be  reconstructed  knowing 
elevations  for  the  1960  Chilean  tsunami.   No  major  tsunamis  were 
generated  in  the  Aleutian-Alaskan  region  from  1857  through  1946;  hence, 
no  reconstruction  is  necessary  for  tsunamis  from  this  region. 
Fortunately,  tsunamis  from  this  region  since  1946  were  representative 
of  the  entire  region  since  these  tsunamis  were  generated  in  the 
western  Aleutians  (1957),  central  Aleutians  (1946)  and  eastern 
Alaskan  area  (1964) . 

The  assumption  that  tsunamis  generated  in  a  single  source 
region  (Kamchatka  or  Chile)  approach,  the  Hawaiian  Islands  from 
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nearly  the  same  direction  is  justified  in  the  case  of  Kamchatka  by 
the  small  spatial  extent  of  the  known  generation  area.   However,  the 
tsunamigenic  area  of  Chile  is  considerably  larger  than  that  of 
Kamchatka.   Still,  the  direction  oP  approach  of  different  tsunamis 
from  Chile  is  nearly  the  same  because  the  tsunamigenic  region  in 
Chile  subtends  a  relatively  small  angle  with  respect  to  the  islands 
as  a  result  of  the  great  distance  between  Chile  and  the  Hawaiian 
Islands  and  the  axis  of  the  Chile  Trench  maintains  a  constant  orienta- 
tion relative  to  the  islands  (thus  directional  radiation  effects  are 
not  important) . 

The  assumption  that  tsunamis  generated  in  a  single  source  region 
have  energy  lying  within  the  same  band  of  wave  periods  is  supported 
by  the  elevation  patterns  produced  along  coastlines  in  the  Hawaiian 
Islands  by  different  historical  tsunamis.   These  patterns  are  such 
that  it  has  been  concluded  that  "tsunamis  of  diverse  geographic 
origin  are  strikingly  different,  whereas  those  from  nearly  the  same 
origin  are  remarkably  similar"  (Eaton,  et  al.,  1961).   A  recent 
study  by  IVybro  (1977)  shows  that  even  the  distributions  of  normalized 
elevations  (i.e.  the  elevation  patterns)  produced  in  the  Hawaiian 
Islands  by  different  Aleutian-Alaskan  tsunamis  are  nearly  the  same 
yet  quite  different  from  the  distributions  for  tsunamis  of  other- 
origins.   This  agreement  occurs  dispite  the  fact  that  the  relatively 
close  proximity  to  Hawaii  of  the  Aleutian-Alaskan  Trench  and  the 
varying  orientation  of  the  trench  axis  relative  to  Hawaii  introduces 
important  directional  effects  for  tsunamis  generated  in  the  Aleutian- 
Alaskan  area.   Apparently,  these  directional  effects  influence  the 


magnitude  of  the  elevations  in  Hawaii  but  do  not  greatly  alter  ele- 
vation patterns.   Therefore,  historical  observations  support  the 
approach  used  in  this  study  which  estimates  the  elevations  produced 
by  tsunamis  from  Chile  or  Kamchatka  prior  to  1946  based  upon  data 
for  tsunamis  from  these  tsunamigenic  regions  recorded  during  the 
years  of  accurate  survey  measurements  from  1946  to  1977. 

In  this  study  elevations  are  calculated  along  the  coastline 
of  the  Hawaiian  Islands  for  the  ten  tsunamis  which  produced  the 
greatest  elevations  in  Hilo,  Hawaii,  from  1S57  through  1977.   The 
finite  element  numerical  model  fills  in  all  data  gaps  for  those  tsunamis 
of  the  ten  that  occurred  from  1946  through  1977  and  then  the  Hilo 
data  (Cox,  1964)  is  used  to  determine  relative  magnitudes  to  reconstruct 
the  elevations  for  those  tsunamis  that  occurred  prior  to  1946.   For 
much  of  the  coast  of  the  Hawaiian  Islands,  the  ten  largest  tsunamis 
since  1837  would  be  the  same  as  the  ten  largest  tsunamis  in  Hilo  for 
the  same  period.   That  is,  the  ten  largest  wouJd  be  the  tsunamis  of 
1960,  1946,  1925,  lS57,'lS77,  1841,  1957,  1952,  and  the  two  in  1868. 
Of  course,  the  order  of  the  ten  largest  would  vary  from  location  to 
location.   For  example,  the  largest  tsunami  on  the  South  Puna  coast 
of  the  big  island  of  Hawaii  was  the  locally  generated  tsunami  of 
1868. 

There  are  a  few  locations  where  historical  data  show  that  one 
of  the  ten  largest  tsunamis  was  not  among  the  ten  largest  in  Hilo 
from  1837  through  1977.   For  example,  the  1896  tsunami  generated 
near  Japan  produced  small  elevations  in  the  Hawaiian  Islands  except 
at  Keauhou  on  the  Kona  coast  of  the  bis  island  of  Hawaii.   The 
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observation  at  Keauhou  was,  therefore,  included  in  the  frequency  of 
occurrence  analysis  for  that  location.   Similarly,  all  other  histor- 
ical observations  of  significant  tsunami  elevations  at  isolated 
locations  were  included  in  the  frequency  of  occurrence  analyses  for 
those  locations,.   Historical  data  for  tsunamis  prior  to  1946  were 
taken  from  Pararas-Carayannis  (1977) .   Data  for  the  local  tsunami  of 
.November  29,  1975,  were  taken  from  the  compilation  by  Loomis  (1976). 
Since  the  local  tsunami  of  1868  had  an  earthquake  epicenter  very 
near  that  of  the  1975  tsunami  and  the  two  apparently  produced  very 
similar  elevation  patterns  (Tilling  et  al . ,  1976),  historical  ele- 
vations for  the  two  tsunamis  recorded  at  the  same  location  were  used  to 
determine  a  relative  magnitude  of  the  tsunamis.  The  many  observations 
for  the  1975  tsunami  were  then  used  to  reconstruct  elevations  for  the 
186S  tsunami. 
Frequency  of  Occurrence  Distribution 

After  the  ten  largest  tsunamis  from  1837  through  1977  were 
determined  at  locations  all  along  the  coasts  of  the  Hawaiian  Islands 
using  the  methods  discussed  in  the  previous  sections,  elevation 
versus  frequency  of  occurrence  curves  were  determined  at  each  location 
by  least-squares  fitting  of  the  data  by  the  following  expression: 

h  =  -B  -  A  log   f  (4.1) 

where 

h  =  elevation  of  maximum  wave  at  the  shoreline, 

f  =  frequency  per  year  of  occurrence. 
Equation  (4.1)  was  used  as  the  frequency  of  occurrence  distribution 
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since  it  has  been  found  to  agree  with  historical  data  at  several 
locations.   Cox  (1964)  found  that  the  logarithm  of  tsunami  frequency 
of  occurrence  was  linearly  related  to  tsunami  elevations  for  the  ten 
largest  tsunamis  occurring  from  1837  through  1964  at  Hilo,  Hawaii. 
Soloviev  (1969]  has  shown  a  similar  relationship  between  tsunami 
frequency  of  generation  and  intensity  for  moderate  to  large  tsunamis. 
Also,  earthquake  frequency  of  occurrence  and  magnitude  have  been 
similarly  related  by  Gutenberg  and  Richter  (1965) ;  IViegel  (1965) 
found  the  same  relationship  for  historical  tsunamis  at  San  Francisco 
and  Crescent  City,  California  and  Adams  (1970)  for  tsunamis  at  Kahuku 
Point,  Oahu.   A  recent  study  by  Rascon  and  Villarreal  (1975)  revealed 
the  same  relationship  for  historical  tsunamis  on  the  west  coast  of 
America,  excluding  Mexico. 

It  is  possible  that  other  distributions  may  agree  with  the 
historical  data  equally  well  as  the  logarithmic  distribution.   For 
example,  the  Gumbel  distribution  has  been  used  in  the  past  to  study 
annual  streamflow  extremes  (Gumbel,  1955).   Borgman  and  Resio  (1977) 
illustrate  the  use  of  this  distribution  to  determine  frequency  curves 
for  non-annual  events  in  wave  climatology.   To  investigate  the  sensi- 
tivity of  calculations  of  100-year  elevations  on  the  assumed  frequency 
distribution,  the  approach  of  Borgman  and  Resio  was  applied  to  the 
Ililo  data  of  Cox.   The  Gumbel  distribution  yields  a  100-year  elevation 
of  28.8  feet  for  data  from  1837  through  1977  and  an  elevation  of  42.5 
feet  for  data  from  1946  through  1977.   This  compares  with  the  elevations 
of  27.3  feet  and  44.2  feet  calculated  for  the  same  time  periods  using 
the  logarithmic  distribution  (Equation  (4.1)).   Clearly,  the  arguments 
used  earlier  concerning  the  period  of  time  that  must  be  considered  for 
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a  valid  analysis  are  not  dependent  upon  the  assumed  frequency  distri- 
bution.  Since  ample  precedent  exists  for  using  the  logarithmic  distri- 
bution for  analysing  tsunami  frequency  of  occurrence,  this  distribution 
is  used  in  this  study. 
Results 

The  methods  described  in  the  previous  sections  allow  frequency  of 
occurrence  curves  to  be  constructed  for  points  all  along  the  coast- 
line of  the  Hawaiian  Islands.   The  coefficients  A  and  B  of  Equation 
(4.1)  contain  all  necessary  information  to  determine  frequency  of 
occurrence  of  tsunami  elevations  at  a  location.   Plots  of  the  coeffi- 
cients A  and  B  versus  location  along  the  coasts  were  constructed  for 
all  of  the  islands.   Smooth  curves  were  drawn  through  these  coeffi- 
cients to  allow  elevation  predictions  at  any  location.   Figures  18 
and  19  show  typical  plots  of  A  and  B  for  the  southeastern  coast  of 
the  big  island  of  Hawaii  (Figure  20).   For  example,  at  Hoopuloa 
(location  11  in  Figure  20)  A  =  6  and  B  =  4;  hence,  the  109-year 
elevation  (F  =  .01)  is  only  8  feet.  The  100-year  elevations  for 
most  locations  in  the  Hawaiian  Islands  are  considerably  larger. 

For  such  situations  involving  human  life,  one  is  less  interested 
in  the  mean  exceedar.ee  frequency  f   (the  average  frequency  per  year 
of  tsunamis  of  equal  or  greater  elevations)  than  in  the  chance  of  a 
given  elevation  being  exceeded  in  a  certain  period  of  time.   To 
calculate  risk,  one  first  notes  that  tsunamis  are  usually  caused  by 
earthquakes,  and  earthquakes  are  often  idealized  as  a  generalized 
Poisson  process  (Ncwmark  and  Rosenblueth,  1971  and  Der  Kiureghian 
and  Ang,  1977).   It  was  assumed  by  Wiegel  (1965)  and  Rascon  and 
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Figure  20.   Location  Map  for  Hawaii. 
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Villarreal  (1975)  that  tsunamis  also  follow  such  a  stochastic  process. 
The  probability  that  a  tsunami  with  an  average  frequency  of  occur- 
rence of  f  is  exceeded  in  D  years,  assuming  that  tsunamis  follow 
a  Poisson  process,  is  given  by  the  following  equation: 

P  =  1  -  e"fD  (4.2) 

Thus,  if  an  acceptable  risk  is  one  chance  in  one  thousand  that  an 
elevation  be  exceeded  during  a  100-year  period,   P  =  .001  and 
D  =  100  years.   Substituting  these  values  in  Equation  (4.2)  yields 
the  elevation  h  which  has  one  chance  in  one  thousand  of  being 
exceeded  during  the  next  100  years. 

West  Coast  of  the  United  States  Predictions 

Introduction 

Unlike  the  Hawaiian  Islands,  the  west  coast  of  the  continental 
United  States  lacks  sufficient  data  to  allow  tsunami  elevation 
predictions  based  upon  local  historical  records  of  tsunami  activity. 
Virtually  all  of  the  west  coast  is  completely  without  data  of  tsunami 
occurrence,  even  for  the  prominent  1964  tsunami.   Only  a  handful  of 
locations  have  historical  data  for  tsunamis  other  than  the  1964  tsunami 
However,  the  Federal  Insurance  Administration  requires  information  on 
tsunami  elevations  for  the  entire  west  coast  of  the  continental  United 
States;  even  for  the  many  locations  that  have  no  known  historical 
data  of  tsunami  activity  and  for  coastal  areas  that  are  currently 
not  developed  (since  these  areas  may  be  developed  in  the  future) . 
The  lack  of  historical  data  of  tsunami  activity  on  the  west  coast 
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necessitates  the  use  of  numerical  models  to  predict  runup  eleva- 
tions. The  Aleutian-Alaskan  area  and  the  west  coast  of  South  .America 
were  found  by  Houston  and  Garcia  (1974)  to  be  the  tsunamigenic 
regions  of  concern  to  the  west  coast  of  the  United  States.   Both 
regions  have  sufficient  data  on  the  generation  of  major  tsunamis  to 
allow  a  statistical  investigation  of  tsunami  generation.   The  genera- 
tion and  deep-ocean  propagation  model  described  in  Chapter  III  can 
then  be  used  to  generate  representative  tsunamis  and  propagate  them 
across  the  deep  ocean.   The  nearshore  numerical  model  also  described 
in  Chapter  III,  can  then  be  used  to  propagate  tsunamis  from  the  deep 
ocean  over  the  continental  slope  and  shelf  to  shore. 

In  this  study,  only  tsunamis  of  distant  origin  are  considered 
in  the  analysis.   Hammack  (1972)  has  shown  that  near  the  generation 
region  of  an  impulsively  generated  tsunami,  the  waveform  is  dependent 
upon  details  of  the  time-dependent  movement  of  the  ground  during  the 
earthquake.   Little  is  known  about  the  actual  time- dependent  ground 
motion  during  earthquakes  generating  tsunamis  and  this  motion  cannot 
be  predicted  in  advance.   Also,  there  is  not  enough  historical  data 
concerning  locally  generated  tsunamis  on  the  west  coast  to  allow 
predictions  of  tsunami  occurrence.   Thus,  reasonable  predictions  on 
the  west  coast  of  the  properties  of  locally  generated  tsunamis  or 
their  likelihood  of  occurrence  are  not  possible  at  this  time. 

The  probability  is  not  considered  very  great  that  a  destructive, 
locally  generated  tsunami  will  occur  on  the  west  coast  of  the  conti- 
nental United  States.   Tsunamis  are  generally  produced  by  earthquakes 
having  fault  movements  that  exhibit  a  pronounced  "dip-slip,"  or 


77 


vertical  component  of  motion.  "Strike-slip/'  or  horizontal  displace- 
ment, fault  movements  are  inefficient  generators  of  tsunamis. 
Faults  on  the  west  coast  of  the  United  States  characteristically 
exhibit  "strike-slip"  motion  since  the  Pacific  block  of  the  earth's 
crust  is  moving, horizontally  relative  to  the  North  American  block. 
The  west  coast  of  the  United  States  does  not  share  the  characteristics 
(ocean  trenches  and  island  arcs)  of  known-generating  areas  and,  in 
fact,  has  not  historically  been  one.   Relatively  small  locally 
generated  tsunamis  have  been  known  to  occur  on  the  west  coast,  but 
there  are  no  reliable  reports  of  major  locally  generated  tsunamis. 
Heights  of  6  feet  in  the  immediate  vicinity  of  the  1927  Point 
Arguello  earthquake  are  the  largest  authenticated  heights  produced 
by  local  tsunamis  on  the  west  coast) .   There  could  be  a  few  locations 
on  the  west  coast  for  which  locally  generated  tsunamis  pose  a  greater 
hazard  than  do  distantly  generated  tsunamis  because  the  elevations 
produced  by  distantly  generated  tsunamis  are  small.   However,  pre- 
dictions of  elevations  produced  by  locally  generated  tsunamis  are 
beyond  the  scope  of  this  study. 
Tsunami  Occurrence  Probabilities 

Historical  data  of  tsunami  generation  must  be  the  basis  for  an 
analysis  that  considers  the  probability  of  tsunami  generation  in  the 
two  tsunamigenic  regions  of  the  Pacific  Ocean  of  concern  to  the  west 
coast  of  the  continental  United  States--the  Aleutian  and  Peru-Chile 
Trench  regions.   A  satisfactory  correlation  between  earthquake 
magnitude  and  tsunami  intensity  has  never  been  demonstrated.   Not 
all  large  earthquakes  occurring  in  the  ocean  even  generate  noticeable 
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tsunamis.   Furthermore,  earthquake  parameters  of  importance  to 
tsunami  generation,  such  as  focal  depth,  rise  time,  and  vertical 
ground  displacement,  have  only  been  measured  for  earthquakes  occurring 
in  recent  years.   Therefore,  data  of  earthquake  occurrence  cannot  be 
used  to  determine  occurrence  probabilities  of  tsunamis.   Instead, 
historical  data  of  tsunami  occurrence  in  generation  regions  must  be 
used  to  determine  these  probabilities. 

In  South  America,  a  wealth  of  information  exists  concerning 
tsunami  generation.   Reliable  data  (grouped  in  intensity  increments 
of  one-half)  exist  for  tsunamis  with  intensity  greater  than  or  equal 
to  0  for  a  171 -year  period  and  greater  than  or  equal  to  2-1/2  for  a 
416-year  period  (Soloviev  and  Go,  1969).   The  intensity  scale  used 
is  a  modification  (by  Soloviev  and  Go)  of  the  standard  Imamura-Iida 
tsunami  intensity.  Intensity  is  defined  as 

I  =  log  (/2  H   )  (4.3) 

This  definition  in  terns  of  an  average  runup  (in  meters)  over  a 
coast  instead  of  a  maximum  runup  elevation  at  a  single  location 
(used  for  the  standard  Imamura-Iida  scale)  tends  to  eliminate  any 
spurious  intensity  magnitudes  caused  by  often  observed  anomalous 
responses  (due,  for  example,  to  local  resonances)  of  single  isolated 
locations. 

Using  the  most  recent  and  complete  catalog  of  tsunami  occurrence, 
in  the  Pacific  Ocean  (Soloviev  and  Go,  1969)  a  relationship  between 
tsunami  intensity  and  frequency  of  occurrence  was  determined  for  the 
tsunamigenic  trench  running  the  length  of  the  Peru-Chile  coast. 
Tsunamis  with  intensity  greater  than  or  equal  to  0  were  considered. 
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It  was  assumed  that  the  logarithm  of  the  tsunami  frequency  of  occur- 
rence was  linearly  related  to  the  tsunami  intensity.   Earlier  in 
this  chapter,  it  was  shown  that  such  a  relationship  holds  between 
both  earthquake  magnitude  and  frequency  of  occurrence  and  tsunami 
intensity  and  frequency  of  occurrence. 

Letting  n(I)   equal  the  probability  of  a  tsunami  with  an 
intensity  1  being  generated  during  any  given  year  and  using 
statistics  for  the  entire  trench  along  the  Peru-Chile  coast,  a 
least-squares  analysis  results  in  the  following  expression: 

n(I)  =  0.74e~°-6iI  (4.4) 

In  using  statistics  for  the  entire  trench  area  along  the  Peru-Chile 
coast,  it  is  assumed  that  the  probability  of  tsunami  occurrence  is 
uniform  along  the  trench.   This  is  a  standard  assumption  for  earth- 
quake frequency  analysis  (Gutenberg  and  Richter,  1965) .   The  tectonic 
justification  of  this  assumption  lies  in  the  fact  that  a  single  Sialic 
block  or  plate  of  the  earth's  crust  or  lithosphere  is  dipping  into 
the  Peru-Chile  Trench  (Wilson,  1959) .   It  can  reasonably  be  expected 
that  the  movement  of  this  single  plate  is  similar  along  its  entire 
length. 

In  the  Aleutian  Trench  region,  only  large  tsunamis  occurring 
in  relatively  recent  years  (since  17SS)  have  been  recorded  due  to  the 
isolation  of  the  area.   Assuming  an  exponential  coefficient  of  -0.71 
for  this  trench  area  (determined  by  Soloviev,  1969)  as  a  mean  value 
for  regions  of  the  Pacific  with  the  most  data  on  tsunamis)  and  using 
only  the  reliable  data  for  large  tsunamis  (intensity  greater  than  or 
equal  to  3.5)  from  Soloviev  and  Go  (1969),  the  following  relationship 
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is  determined  by  a  least-squares  analysis: 

nl    n  it-  -0.711 

n  ( I  j  =  0 . 1 1 3  e 
Again,  the  probability  of  tsunami  occurrence  is  assumed  to  be  uniform 
along  the  trench. 
Use  of  Deterministic  Numerical  Models 

To  relate  the  probability  distributions  of  tsunami  intensities  to 
source  characteristics,  it  is  assumed  that  the  ratio  of  the  source 
uplift  heights  producing  two  tsunamis  of  different  intensity  (as 
defined  in  the  previous  section)  is  equal  to  the  ratio  of  the  average 

runup  heights  produced  on  the  coasts  near  these  tsunami  sources.  This 

|lrI2l 
ratio  is  equal  to  2         for  two  tsunamis  with  intensities   I  and 

l2 

If  H  is  the  wave  height  in  the  direction  parallel  to  the  maior 
a 

axis  of  length  a  cf  a  tsunami  source  with  an  elliptical  shape 

(large  tsunamis  have  historically  had  elliptically-shaped  uplifts) 

and  H.   is  the  wave  height  in  the  direction  parallel  to  the  minor 
b 

axis  of  length  b,   then  experimental  research  of  tsunami  generation 

has  shown  that  H, /H   is  approximately  equal  to  a/b   (Hatori, 

b  a 

1963).  For  a  large  tsunami.  H,   can  be  larger  than  H   by  a  factor 

J  a  k  o  g 

of  as  much  as  5  or  6.   Thus,  the  orientation  of  the  tsunami  source 
relative  to  the  area  where  runup  is  to  be  determined  is  very  important, 
i.e.,  the  runup  at  a  distant  site  due  to  the  generation  of  a  tsunami 
at  one  location  along  a  trench  cannot  be  considered  as  being  repre- 
sentative of  all  possible  placements  of  the  tsunami  source  in  the 
entire  trench  region.   Hence,  the  Aleutian  and  Peru-Chile  Trenches 
had  to  be  segmented  and  runup  along  the  west  coast  of  the  United 
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States  determined  for  tsunami  sources  located  at  the  center  of  each 
of  the  segments. 

The  spatial  size  of  a  tsunami  source  was  standardized  since 
there  is  not  an  apparent  correlation  between  tsunami  intensity  and 
spatial  size  of  a  tsunami  source.   for  example,  the  1946  Aleutian 
tsunami  had  an  uplift  region  of  very  small  spatial  extent,  whereas 
the  1957  Aleutian  tsunami  had  an  uplift  region  that  covered  perhaps 
the  greatest  spatial  extent  of  any  known  earthquake  (Kclleher, 
1972) .   Yet  the  1946  tsunami  had  the  greater  intensity,  producing, 
in  general,  greater  runup  elevations  in  the  near  and  distant  regions. 
The  standard  source  used  is  discussed  in  detail  by  Houston  and 
Garcia  (1974).   It  is  elliptical  in  plan  view  with  a  major  axis  length 
of  600  miles  and  a  minor  axis  length  of  130  miles.   The  vertical  dis- 
placement increases  linearly  from  a  zero  elevation  on  the  ellipse 
perimeter  to  a  maximum  displacement  of  30  feet.   The  vertical  dis- 
placement has  a  parabolic  crest  (concave  downward)  parallel  to  the 
direction  of  the  major  axis.   The  standard  source  represents  a  large 
tsunami  with  intensity  4  on  the  modified  Imamura-Iida  scale. 
Certainly,  tsunamis  of  low  intensity  may  have  smaller  spatial  extents; 
however,  large  tsunamis  pose  the  greatest  threat  to  a  distant  area 
such  as  the  west  coast  of  the  United  States.   These  large  tsunamis 
can  be  expected  to  have  similar  spatial  extents,  with  any  spatial 
differences  being  unimportant  in  the  far- field  compared  with  the  effects 
of  source  orientation  and  vertical  uplift.   Vertical  uplifts  are  assigned 
to  different  intensities  in  accordance  with  the  convention  discussed 
at  the  beginning  of  this  section. 


82 


Figures  21  and  22  show  the  Aleutian  Trench  divided  into  12  seg- 
ments and  the  Peru-Chile  Trench  into  3  segments.   The  segments  in  the 
Aleutian  Trench  were  approximately  one-quarter  the  length  of  the  major 
axis  of  the  standard  source,  whereas  the  segments  in  the  Peru-Chile 
Trench  were  approximately  the  length  of  the  major  axis  of  the  standard 
source.   The  standard  source  was  centered  in  each  segment  such  that 
the  major  axis  of  the  source  was  parallel  to  the  trench  axis.  Uplift 
regions  historically  have  had  such  an  orientation  relative  to  trench 
systems.   The  Aleutian  Trench  is  segmented  much  finer  than  the  Peru- 
Chile  Trench  since  the  Aleutian  Trench  is  oriented  relative  to  the 
west  coast  such  that  elevations  produced  on  the  west  coast  are  very 
sensitive  to  the  exact  location  of  a  source  along  the  Trench.   Uplifts 
along  the  Peru-Chile  Trench  do  not  radiate  energy  directly  toward  the 
west  coast  regardless  of  their  position  along  the  Trench.   The  Peru  and 
Chile  sections  of  the  Peru-Chile  Trench  have  constant  orientations 
relative  to  the  west  coast  of  the  United  States;  therefore,  eleva- 
tions on  the  west  coast -are  not  very  sensitive  to  source  location 
within  these  sections. 

In  each  of  the  segments  of  the  Aleutian  and  Peru-Chile  Trenches 
tsunamis  with  intensities  from  2  to  5  in  steps  of  one-half  intensity 
are  generated  and  propagated  across  the  deep  ocean  using  the  deep- 
ocean  numerical  model  discussed  in  Chapter  II.   Tsunamis  with  inten- 
sity less  than  2  are  too  small  to  produce  significant  runup  on  the 
west  coast.   An  upper  limit  of  5  was  chosen  because  the  largest 
tsunami  intensity  ever  reported  was  less  than  5  (Soloviev  and  Go, 
1969).   Gutenberg  and  Richter  (1965)  indicate  that  there  is  an  upper 
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limit  to  the  strain  that  can  be  supported  by  rock  before  fracture. 
Thus,  earthquakes  only  reach   certain  maximum  magnitudes  and  tsunamis 
can  be  expected  to  have  similar  upp . r  limits  to  intensity.   Perkins 
(1972)  and  McGarr  (1976)  have  demon ^tro Led  that  future  earthquakes 
cannot  have  seismic  moments  (measure  of  earthquake  magnitude  for 
large  earthquakes)  much  larger  than,  those  of  earthquakes  that  occurred 
in  recorded  history. 

The  waveforms  propagated  to  the  west  coast  by  the  deep-ocean 
model  are  used  as  before  as  input  to  the  nearshore  propagation  model. 
Each  waveform  is  propagated  from  a  water  depth  of  500  meters  to 
shore  using  the  nearshore  model  and.  one  of  the  grids  shown  in  Figure 
13.   Thus,  at  each  grid  location  on  the  shoreline  of  the  west  coast, 
there  is  a  group  of  105  waveforms  -  seven  waveforms  (for  intensities 
from  2  to  5  in  one- half  intensity  increments)  for  each  segment,  of 
the  Aleutian  and  Peru-Chile  Trenches.  Each  of  these  waveforms  has  an 
associated  probability  equal  to  the  probability  that  a  certain 
intensity  tsunami  will  be  generated  in  a  particular  segment  of  a 
trench  region. 
Effect  of  the  Astronomic al  Tides 

The  maximum  "still-water"  elevation  produced  during  tsunami 
activity  is  the  result  of  a  superposition  of  tsunamis  and  tides. 
Therefore,  she  statistical  affect  of  the  astronomical  tides  on 
total  tsunami  runup  must  be  include.?  in  the  predictive  scheme  pre- 
sented in  this  study.   Since  the  waveforms  calculated  by  the  nearshore 
model  go  not  have  a  simple  form  (e.g.  sinusoidal),  the  statistical 
effect  of  the  astronomical  tide  on  tsunami  runup  must  be  determined 
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through  a  numerical  approach.   (The  effect  of  the  astronomical  tides 
was  not  considered  in  the  section  of  predictions  for  the  Hawaiian 
Islands  because  the  tidal  range  is  quite  snail  for  these  islands  and 
the  local  historical  data  implicitly  contain  the  effects  of  the 
tides. ) 

The  waveforms  calculated  by  the  nearshore  numerical  model 
extend  over  a  period  of  time  of  approximately  two  hours.   Three  or 
four  wave  crests  (the  largest  waves  in  the  tsunami)  arrive  during 
this  time.  Smaller  waves  arriving  at  later  times,  however,  often 
persist  for  days  at  a  coastal  location.  An  analysis  of  tide  gage 
records  of  the  1960  and  1964  tsunamis  on  the  west  coast  indicates 
that  these  smaller  waves  have  amplitudes  on  the  average  of  40  percent 
of  the  maximum  wave  amplitude  of  the  tsunami.   Therefore,  a  sinu- 
soidal group  of  these  smaller  waves  were  added  to  each  of  the  calculated 
waveforms  so  that  the  total  waveform  extended  over  a  24  hour  period. 
These  smaller  waves  are  important  for  locations  where  tsunamis 
are  fairly  small  compared  with  tidal  variations.   At  such  locations 
the  maximum  combined  tsunami  and  astronomical  elevation  occurs 
during  the  maximum  tidal  elevation. 

A  computer  program  was  developed  to  predict  time-histories  of 
the  astronomical  tides  at  all  grid  locations  on  the  west  coast  of 
the  United  States.   The  program  was  based  upon  the  harmonic  analysis 
methods  used  in  the  past  by  the  Coast  and  Geodetic  Survey  for  mechan- 
ical tide-predicting  machines  (Schureman,  1948).   Tidal  constants 
available  from  the  Coast  and  Geodetic  Survey  were  used  as  input  to 
the  computer  program;  and  a  year  of  tidal  elevations  was  then  predicted 
for  grid  locations  all  along  the  west  coast.   The  year  1964  was 
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selected  because  all  the  major  tidal  components  had  a  node  factor 
equal  to  approximately  1.00  during  this  year  making  it  an  average 
year  (node  factors  can  vary  by  about  10  percent).   The  node  factor 
is  associated  with  the  revolution  of  the  moon's  node  and  has  an 
18.6-year  cycle.   Since  a  tsunami  can  arrive  at  any  time  during  this 
18.6-year  period  (arrival  at  a  low  of  the  node  factor  is  equally 
likely  as  an  arrival  at  a  high],  the  statistical  effect  of  the 
varying  node  factor  is  small  and  an  average  value  should  be  used. 
The  statistical  effect  of  the  varying  node  factor  on  the  predicted 
runup  elevations  can  be  shown  to  be  a  small  fraction  of  an  inch, 
using  the  approach  discussed  in  the  next  section  and  Appendix  A 
(with  the  variance  of  the  nodal  variation  equal  to  approximately 
0.1  feet  squared). 

The  year  of  tidal  elevations  calculated  at  each  of  the  nearshore 
numerical  model  grid  points  along  the  west  coast  was  then  subdivided 
into  15-minute  segments.   The  24-hour  waveforms  were  allowed  to 
arrive  at  the  beginning  of  each  of  these  15-minute  segments  and 
then  superposed  upon  the  astronomical  tide  for  the  24-hour  period. 
The  maximum  combined  tsunami  and  astronomical  tide  elevation  over 
the  24-hour  period  was  determined  for  tsunamis  arriving  during  each 
of  these  15-minute  starting  times  during  a  year.   All  of  the  maximum 
elevations  had  an  associated  probability  equal  to  the  probability 
that  a  certain  intensity  tsunami  would  be  generated  in  a  particular 
segment  of  the  two  trench  regions  and  arrive  during  a  particular  15- 
minute  period  of  a  year. 


The  many  maximum  elevations  with  associated  probabilities  can 
be  used  to  determine  cumulative  probability  distributions  of  combined 
tsunami  and  astronomical  tide  elevations.   The  maximum  elevations 
are  ordered  and  probabilities  summed,  starting  with  the  largest 
elevations,  until  a  desired  probability  is  obtained.   The  elevation 
encountered  when  the  summed  probabilities  reach  a  desired  value  P 
is  the  elevation  that  is  equaled  or  exceeded  with  an  average  fre- 
quency of  once  every  1/P  years.   Thus,  when  the  summed  probabilities 
reach  the  value  .01,  the  elevation  associated  with  the  last  prob- 
ability summed  is  the  100-year  elevation. 

The  100-year  and  500-year  elevations  are  determined  at  all 
grid  points  of  the  nearshore  numerical  grids  using  techniques  des- 
cribed previously  in  this  chapter.   Smooth  curves  are  then  used  to 
connect  all  discrete  elevations  so  that  continuous  predictions  of 
elevation  versus  coastline  location  can  be  made  for  all  the  west 
coast.   Figure  25  shows  100-year  and  500-year  elevations  in  the 
Crescent  City,  California,  area. 
Comparison  with  Local  Observation  Predictions 

Crescent  City  and  San  Francisco,  California,  are  the  only 
locations  on  the  portion  of  the  west  coast  of  the  United  States 
considered  in  this  study  that  have  sufficient  historical  data  of 
tsunami  activity  to  allow  frequency  of  occurrence  predictions  based  upon 
local  historical  observations.   Wiegel  (1965)  made  such  predictions  of 
tsunami  height  (trough  to  crest  height)  for  Crescent  City  based  upon 
the  period  from  1900  to  1965.   He  predicted  100-year  and  500-year 
heights  at  Crescent  City  of  approximately  25.6  feet  and  43.2  feet, 
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respectively.   If  the  crest  amplitude  is  taken  to  be  one-half  the  total 
height,  the  100-year  and  500-year  elevations  are  12.8  feet  and  21.6 
feet,  respectively.   However,  the  crest  amplitude  at  Crescent  City  is 
typically  greater  than  one-half  the  wave  height  (e.g.  the  crest  ampli- 
tude of  the  largest  wave  of  the  1964  tsunami  was  approximately  60  per- 
cent of  the  total  height).   If  Wiegel's  analysis  is  applied  to  historical 
crest  elevations  instead  of  heights,  the  100-year  and  500-year  elevations 
are  found  to  be  15.4  and  26.4  feet,  respectively.   Furthermore,  the 
analysis  can  now  be  applied  to  the  longer  time  period  from  1900  through 
1977.   The  100-year  and  500-year  elevations  based  upon  this  longer  time 
span  are  14.5  feet  and  25.5  feet  respectively. 

Figure  23  shows  100-year  and  500-year  elevations  of  13.1  feet 
and  24.9  feet,  respectively.   These  values  compare  very  favorably 
with  the  14.5  feet  and  25.5  feet  elevations  determined  from  historical 
data  for  the  period  of  time  from  1909  through  1977.   The  elevations 
predicted  by  the  analysis  based  upon  the  local  historical  data  are 
probably  somewhat  larger  than  the  elevations  predicted  in  this 
report  because  the  short  time  period  (relative  to  100  or  500  years) 
from  1900  through  1977  includes  the  exceptionally  active  years  of 
tsunami  generation  from  1946  through  1964  (see  Chapter  III) . 

The  13.1-foot  and  24.9-foot  elevations  predicted  for  Crescent 
City  using  the  techniques  discussed  in  this  study  are  not  totally 
comparable  to  the  elevations  predicted  using  Wiegel's  analysis  based 
upon  local  historical  data  since  the  effect  of  the  astronomical  tide 
has  been  included  in  the  elevations  predicted  in  this  study.   However, 
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the  statistical  effect  of  the  astronomical  tide  on  the  total  ele- 
vation is  not  significant  at  Crescent  City  due  to  the  large  amplitude 
of  tsunamis  there.   Appendix  A  shows  that  the  statistical  effect  of 

the  astronomical  tide  for  a  location  where  tsunamis  are  large  is  to 

2 
ceo 
increase  the  predicted  elevation  by  an  amount  equal  to  — x—  , 

2  °°   2 

where  a   is  the  tidal  variance  and  equals   E  C   ,   C   is 

1  m     m 

m=l 

equal  to  the  m   tidal  constituent,  and  a  is  given  by  the  following 
expression  [since  it  was  shown  earlier  in  this  chapter  that  elevations 
are  linearly  related  to  the  logarithm  of  probabilities): 

P(Z)  -  Ae~a~ 
with 

Z  =  the  elevation  above  local  mean  sea  level,  and 
P(Z)  =  the  cumulative  probability  distribution  for  the  elevation 
at  a  given  site  being  equal  to  or  exceeding   Z  due  only  to  the 

maximum  wave  of  the  tsunami.   Based  upon  the  tidal  constituents 

2 
predicted  by  the  Coast  and  Geodetic  Survey,  o  =7.1  fox  Crescent 

City.   Using  Wiegel's  data  for  the  period  of  time  from  1900  through 

1977  gives  a  =  .145   .   Therefore,  the  astronomical  tide  contributes 

approximately  0.5  feet  to  the  elevations  predicted  for  Crescent 

City. 

Frequency  of  occurrence  calculations  just  inside  San  Francisco 

Bay  were  made  in  an  earlier  study  (Houston  and  Garcia,  1975)  using  a 

simple  one- dimensional  analytical  solution  for  nearshore  propagation. 

However,  portions  of  San  Francisco  Bay  are  included  in  the  nearshore 

numerical  model  (see  Figure  13)  used  in  this  study  in  order  that  the 

effect  of  the  Bay  on  elevations  outside  the  Bay  can  be  properly 
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simulated.  Thus,  approximate  elevation  calculations  can  be  made  in 
this  study  at  the  location  of  the  Presidio  tide  gage  within  San 
Francisco  Bay.  Figure  24  shows  the  100-year  and  500-year  elevation 
predictions  at  the  Presidio  to  be  6.1  feet  and  10.0  feet,  respectively. 

Local  historical  data  of  tsunami  activity  at  the  Presidio  was 
compiled  by  Wiegel  (1965).   Using  a  logarithmic  distribution  and 
Wiegel's  wave  height  data  for  San  Francisco  over  the  period  1900 
through  1977  yields  100-year  and  500-year  heights  (trough  to  crest) 
of  8.0  feet  and  14.2  feet  respectively.   If  the  crest  amplitude  is 
one-half  the  height,  this  yields  100-year  and  500-year  elevations  of 
4.0  feet  and  7.1  feet  respectively.   The  same  analysis  based  upon 
data  of  crest  amplitudes  yields  100-year  and  500-year  elevations  of 
4.4  feet  and  7.8  feet,  respectively. 

It  is  difficult  to  compare  the  elevations  predicted  using  local 
historical  data  with  the  elevations  predicted  in  this  study  for  San 
Francisco,  since  the  effect  of  the  astronomical  tides  (included  in 
the  calculations  of  this  study)  cannot  easily  be  estimated,  as  was 
the  case  for  Crescent  City.   Tsunamis  recorded  at  the  Presidio  in 
San  Francisco  are  known  to  persist  at  fairly  substantial  levels  for 
extended  periods  of  time  apparently  due  to  an  oscillation  phenomenon. 
For  example,  during  the  1964  tsunami  there  were  six  waves  of  approxi- 
mately equal  or  greater  amplitude  than  the  initial  wave  (in  addition 
to  smaller  waves)  during  the  first  8  hours  of  the  tsunami.   Oscilla- 
tions approximately  40  percent  of  the  amplitude  of  the  initial  wave 
persisted  for  at  least  another  24  hours.   The  tidal  range  at  San 
Francisco  is  5.7  feet  with  mean  higher  high  water  (MHHW)  2.7  feet 
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above  mean  sea  level  (MSL)  and  mean  high  water  (MHW)  1.0  feet  above 
MSL.   At  least  one  of  the  main  waves  of  a  tsunami  at  San  Francisco 
probably  superposes  upon  either  mj iiv'  (with,  the  smaller  oscillations 
then  adding  to  MHHW)  or  MJIIIW.   Thus,  the  effect  of  the  tides  is  to 
contribute  1  to. 2. 7  feet  to  the  total  elevation.   These  elevations 
when  added  to  the  elevations  predicted  from  local  historical  data 
yield  100-year  and  500-year  elevation:;  of  5.4  feet  to  7.1  feet  and 
8.8  feet  to  10.5  feet,  respectively.   These  elevations  compare 
favorably  with  the  6.1  feet  and  10.0  feet  elevations  determined  in 
this  study. 


CHAPTER  V 
CONCLUSIONS 

The  major  objective  of  the  present  study  has  been  to  determine 
tsunami  elevation  frequencies  of  occurrence  at  the  shoreline  for  the 
Hawaiian  Islands  and  the  west  coast  of  the  continental  United  States. 
Numerical  models  were  developed  to  simulate  nearshore  tsunami  propa- 
gation and  interaction  with  coastlines.   A  finite  element  numerical 
model  was  used  to  interpolate  between  modern  historical  measurements 
of  tsunami  elevations  in  the  Hawaiian  Islands.   Finite  difference 
numerical  models  were  used  to  propagate  tsunamis  from  the  Aleutian- 
Alaskan  and  Peru-Chile  regions  to  the  west  coast  of  the  United  States 
so  that  frequency  of  occurrence  predictions  could  be  made. 

The  good  agreement  between  tide  gage  recordings  and  numerical 
simulations  of  the  1960  Chilean  and  1964  Alaskan  tsunamis  in  the 
Hawaiian  Islands  indicates  that  linear,  nondispersive,  and  dissipation- 
free  equations  govern  nearshore  propagation  of  the  main  leading  waves  of 
long  period  tsunamis  for  coastal  regions  having  very  short  continental 
shelves.   A  short  continental  shelf  limits  the  time  available  for  the 
development  of  effects  governed  by  the  terms  neglected  in  the  equation 
of  motion.   Long  wave  equations  that  include  nonlinear  and  bottom 
stress  effects  were  shown  to  govern  propagation  of  the  main  leading 
waves  of  long  period  tsunamis  over  the  shallow  water  region  of  the  west 
coast  of  the  United  States. 
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Good  agreement  between  the  numerical  model  simulations  and  tide 
gage  recordings  of  historical  tsunamis  confirms  the  experimental  flume 
tests  by  Hammack  (1972)  that  showed  that  the  far-field  waveform  generated 
by  a  large  impulsive  ground  displacement  does  not  depend  upon  the 
transient  motions  accompanying  the  permanent  ground  motion.   The  good 
agreement  also  shows  that  the  accepted  permanent  ground  deformations  of 
the  1960  Chilean  and  1964  Alaskan  tsunamis,  which  were  partly  measured 
and  partly  estimated  by  Plafker  and  Savage  (1970)  and  Plafker  (1969) 
respectively,  are  fairly  accurate. 

Finally,  this  study  shows  that  numerical  models  can  be  used  in 
conjunction  with  limited  available  historical  information  to  make 
reasonable  predictions  of  tsunami  elevations  on  coastlines.   For  an 
area  having  local  historical  data  of  tsunami  occurrence,  such  as  the 
Hawaiian  Islands,  numerical  models  can  be  used  to  predict  elevations  at 
particular  locations  not  having  data.   For  an  area  lacking  almost  any 
local  historical  data,  such  as  the  west  coast  of  the  United  States, 
elevation  predictions  can  be  made  based  upon  data  of  tsunami  occurrence 
in  tsunamigenic  regions  and  numerical  model  propagation  of  tsunamis  from 
generation  regions  to  the  area  of  interest. 


APPENDIX 

SUPERPOSITION  OF  LARGE  TSUNAMIS  AND  ASTRONOMICAL  TIDES 

If  the  largest  wave  of  a  tsunami  produces  the  largest  total  ele- 
vation when  the  tsunami  is  superposed  upon  the  astronomical  tide 
(regardless  of  where  the  largest  wave  falls  in  the  tidal  cycle),  the 
statistical  effect  of  the  astronomical  tide  on  the  total  elevation 
can  be  estimated  analytically.   Tsunamis  satisfy  this  condition  when 
their  maximum  waves  have  amplitudes  larger  than  the  sum  of  the  maximum 
tidal  elevation  and  the  elevation  of  one  of  the  smaller  waves  which 
follow  the  major  waves  and  persist  over  periods  of  days. 

Let 

Z  =  the  elevation  at  any  time  above  local  mean  sea  level. 

P„(Z)  =  the  cumulative  probability  distribution  for  the 

elevation  at  a  given  site  being  equal  to  or  exceeding 
Z  due  only  to  the  maximum  wave  of  the  tsunami. 

P  (Z)  =  the  probability  of  the  elevation  at  the  same  location 
being  equal  to  or  exceeding  Z  clue  only  to  the 
astronomical  tide. 

P(Z)  =  the  cumulative  probability  distribution  for  the  eleva- 
tion at  a  given  site  being  equal  to  or  exceeding  Z 
due  to  the  maximum  wave  of  the  tsunami  and  the 
astronomical  tide. 

According  to  Chandrasekhar  (1943),  P(Z)  can  be  calculated  from  the 

equation 

P(Z)  =   /   fg(A)Ps(Z-A)dX  (Al) 
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where 


dPR(Z) 


'fr  '  dz 

and  ffi^  is  the  Probabilit>'  density  for  the  astronomical  tide. 

Pg(Z)  can  be  approximated  by  a  Gaussian  distribution. 
Therefore, 

c   cn>  1    -Z  /2a 

fB^)  -  -—  e  (A2) 

/2ir  c 

where  the  variance  is 

2    °°    2 

a  =   E  C 

and  C   equals  the  m   tidal  constituent. 

PS(Z)  can  be  approximated  by  an  exponential  function  (see 
Chapter  II)  as  follows: 

PSC)  =  Ae~uZ  (A3) 

Substituting  Equations  (A2)  and  (A3)  into  (Al)  and  performing  the 
integration  yields 

P(Z)  =  Ae 


2 
■a(Z =-) 


71 
or  p(z)  =  Ae"a"J 

where  z1  =  Z  -  °°-- 

2 

Thus,  the  net  effect  of  the  astronomical  tide  is  to  produce  a  P(Z) 

2 
identical  with  P„(Z)  except  for  a  shift  of  Z  by  an  amount  ^~ 
o  '  2 

2 
The  variance  a     can  be  determined  for  locations  in  the  United 

States  from  tidal  harmonic  coefficients  available  from  the  Coast 

and  Geodetic  Survey.   The  exponential  constant  a   in  Equation  (A3) 

is  obtained  from  a  least-squares  exponential  fit  of  tsunami  elevation 

data  at  a  location. 
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